The antagonist yeast K. apiculata 34-9 was isolated and preserved in a postharvest laboratory of Huazhong Agricultural University. It was cultured on YPDA (1% yeast extract, 2% peptone, 2% dextrose, 2% agar) plate for 2 days, or YPD medium (1% yeast extract, 2% peptone, 2% dextrose) in a rotary shaker coupled with 150 rpm shaking for 1 day at 25°C, respectively.

The plant pathogen P. digitatum N1 which was isolated and preserved in a postharvest lab of Huazhong Agricultural University was used in this study. An even mycelial layer was obtained according to a previous study (Tian et al., 2020). Briefly, 1 mL spore suspensions at a concentration of 1.0 × 10⁶ CFU/mL of P. digitatum was transferred into 100 mL liquid PDA medium (below 50°C) and immediately divided into Petri dishes followed by absolutely mixed. The Petri dishes were air-dried, then incubated in a climate chamber for 4 days at 25°C.

The Hi-C sequencing approach was explored to obtain high-quality assembly and annotated genome. Fresh K. apiculata 34-9 cell samples grown on YPDA plates were harvested after incubation for 2 days at 25°C, then fixed in 1.0% formaldehyde at room temperature for crosslinking to maintain its cellular three-dimensional structure. The yeast cells were lysed, and the cross-linked DNA obtained was then digested with HindIII restriction enzyme. The digested DNA ends were filled and labeled with biotin, then ligated to form a circular chimeric molecule, and the DNA was purified and cut into 300–700 bp fragments. The density and insert size of the constructed library was determined using Qubit 2.0 (Thermo Fisher Scientific Inc., Waltham, United States) and Agilent 2100 Bioanalyzer system (Agilent Technologies Inc., Santa Clara, United States). Then the constructed library was processed for paired-end sequencing on the Illumina Hiseq X Ten platform. The Hi-C reads were filtered using Hi-C Pro to validate paired-end reads (Servant et al., 2015). In addition, the draft contigs obtained from Hi-C sequencing were grouped, ordered, and oriented for chromosomes assembly using the software Lachesis with default parameters (Burton et al., 2013).

Tandem repeats were identified with the utilization of a Tandem Repeat Finder (TRF) in the genome of K. apiculata 34-9 (Benson, 1999). Transposon elements (TEs) of K. apiculata 34-9 genome were annotated with the combination of homology and de novo approaches. Briefly, RepeatMasker¹ and RepeatProteinMask were used for the annotation of TEs of K. apiculata 34-9 genome based on the RepBase database² and RepeatMasker package TE protein database, respectively. RepeatModeler and LTR_FINDER were applied for de novo repeat prediction (Xu and Wang, 2007), and RepeatMasker was utilized for the annotation of TEs with default parameters.

De novo gene prediction and homology-based strategy were employed to perform the gene prediction and functional annotation, and these results were integrated into a gene model using MAKER2 (Holt and Yandell, 2011). The completeness and accuracy of the assembly was assessed using CEGMA (Parra et al., 2007), and the final reliable gene set was obtained in a HiCESAP process. Protein sequences of H. guilliermondii (accession number: GCA_900119595.1), H. opuntiae (accession number: GCA_001749795.1), H. uvarum (accession number: GCA_001747055.1), H. valbyensis (accession number: GCA_001664025.1) genomes were aligned to the genome of K. apiculata 34-9 to perform homology-based gene prediction. All predicted protein-coding sequences were aligned to the SwissProt, non-redundant protein (NR), Gene Ontology (GO), InterPro, Kyoto Encyclopedia of Genes and Genomes (KEGG), and TrEMBL databases for further annotation using BLASTP with default parameters.

For annotation of non-coding RNA, tRNAscan-SE was used to screen the tRNA of K. apiculata 34-9 genome, and BLASTN was applied to identify rRNA based on sequence alignment with homologous species, micro-RNA (miRNA) and small nuclear RNA (snRNA) were annotated by INFERNAL from Rfam.³

The re-assembled genome of K. apiculata 34-9 was compared with H. uvarum (assembly GCA_001747055.1) and S. cerevisiae (assembly GCA_002057635.1). Collinearity among these three yeasts was performed using TBtools (Chen et al., 2020) with default parameters.

For low temperature (LT) treatment, yeast samples were collected after incubation for 1 day in YPD at 25°C, the cells were rapidly frozen in liquid nitrogen after centrifugation at 3,000 × g for 10 min, and marked as 0 h, then stored at −80°C. From that time, some YPD flasks were placed into a 4°C incubator for continuous incubation, and yeast samples were harvested after 5, 12, 24, and 48 h, respectively. Yeast grown at 25°C (normal temperature, Kap) was handled as the same procedure.

For the preparation of interaction between antagonist K. apiculata 34-9 and P. digitatum N1 (Pdi), yeast samples harvested after incubation for 1 day in YPD at 25°C were adjusted to a concentration of 1.0 × 10⁸ CFU/mL, and a total of 1 mL K. apiculata cell suspension was inoculated to each P. digitatum plate and continuously cultured at 25°C, from that time (marked 0 h), samples were harvested as the above time course. There were three biological replicates of each sample. All samples were rapidly frozen in liquid nitrogen and immediately stored at −80°C.

Total RNA from collected yeast samples was isolated using TRIzon Reagent (Cwbiotech, China) following the manufacturer’s instruction. RNA concentration and integrity were measured using NanoDrop 2000 (Thermo Fisher Scientific Inc., Waltham, United States) and Agilent 2100 Bioanalyzer system (Agilent Technologies Inc., Santa Clara, United States), respectively. RNA-sequencing (RNA-seq) libraries were prepared using a NEBNext Ultra RNA Library Prep Kit for Illumina (New England BioLabs^® Inc.) according to the manufacturer’s guide, and the libraries were sequenced on an Illumina HiSeq X 10 platform (Illumina Inc., San Diego, CA). The generated reads were trimmed to eliminate adaptors and enhance quality. The filtrated clean reads were mapped to K. apiculata 34-9 genome using Hisat2 (version 2.1.0) and aligned to the generated transcriptome assembly using Bowtie2 (version 2.2.2), respectively. The output results were estimated using RSEM (version 1.3.0) to obtain normalized counts of gene expression level. The values of expression genes were calculated as fragments per kilobase per million mapped reads (FPKM).

The analyses of differentially expressed genes (DEGs) were conducted using R package DEseq2, and the DEGs were selected if their false discovery rate (FDR) < 0.05 and a log2 fold change > 1, then the DEGs were annotated to GO terms (Ye et al., 2018) and KEGG enrichments (Xie et al., 2011) for functional analyses.

Quantitative real-time PCR (qRT-PCR) assays were conducted to confirm RNA-seq results independently. Total RNA was extracted with Fungal RNA Kit (OMEGA Biotech, Guangzhou, China) followed by the operating instructions. The first-strand cDNA synthesis was conducted using a PrimeScript^TM RT reagent Kit with gDNA Eraser (TaKaRa Biotechnology, Dalian, China). The PowerUP^TM SYBR^TM Green Master Mix (Life Technologies, United States) was used for the detection of the relative expression profiles according to the recommended procedures. qRT-PCR was performed on a Q6 Flex Real-Time PCR system (Life Technologies, United States). The relative gene expression level was analyzed by the 2^–ΔΔCt method (Livak and Schmittgen, 2001). The β-actin was used as a reference gene for normalization, all qRT-PCR primers were designed using Primer-BLAST⁴ and shown in Supplementary Table 1.

The Hi-C sequencing and RNA-seq data generated in this study were released in NCBI Sequence Read Archive database⁵ under accession numbers PRJNA667993 and PRJNA668195.

The next-generation sequencing (previous data) and Hi-C scaffolding were combined to improve K. apiculata yeast genome assembly and chromosome anchoring accuracy. A total of 3.47 Gb (428.40 × coverage) clean data were generated after quality control. K. apiculata genome was assembled into 8.1 Mb of scaffolds sequences and the assembly contained 15 scaffolds (more than 1 kb per scaffold), with a scaffold N50 of 1,349,021 bp (about 4.8-fold larger than that of the previous version), and the maximum scaffold length was 1,727,301 bp (Table 1). Meanwhile, the assembly comprised 155 contigs (more than 1 kb per contig), with a contig N50 of 104,210 bp, and the maximum contig length was 553,666 bp. The GC content was 31.98%, and there were 4,014 protein-coding genes of the new genome version, which was 228 more than the previous annotation version. In addition, 99.51% (8.06 Mb) of the genome sequences had been divided into seven groups (Supplementary Table 2). As many as 75.55% of the reads were uniquely mapped to the genome of K. apiculata 34-9, which were good enough for further analyses. A Hi-C contact heatmap was visualized to confirm the accuracy of Hi-C assembly at 20 kb resolution, all bins could be clearly anchored onto seven pseudochromosomes, the intensity of red color represents the number of Hi-C links on the pseudochromosomes of final assembly (Figures 1A,B). In addition, synteny blocks among K. apiculata 34-9, H uvarum and S. cerevisiae were visualized (Figure 1C).

Repeat elements are important parts of the genome. A total of 133,569 bp (1.65%) repeat sequences were identified in the genome of K. apiculata 34-9. Among them, the proportion of long terminal repeats (LTRs), long interspersed repeats, short interspersed repeats, and DNA transposon was 0.58, 0.16, 0.01, and 0.36% of the genome, respectively (Supplementary Table 3).

In the present study, we aimed to uncover the mystery of K. apiculata responses to cold shock and P. digitatum stress, and comparative transcriptomics studies were conducted to gain more insights into the underlying molecular basis. A total of 304.08 Gb clean data were generated from the sequenced 39 yeast samples (at least 21,559,868 reads pairs of each sample), the Q30 was more than 92.54 and 80.54–85.50% of the reads were uniquely mapped to the genome of K. apiculata 34-9 (Supplementary Table 4). Since the trials were performed within different time courses, yeast samples obtained at the initial time (incubation for 1 day at 25°C) were marked as 0 h, and their transcriptome data were filtered out, then samples from the rest of the four time periods were conducted for transcription profiles analyses.

The detailed number of DEGs are illustrated in Table 2, and as many as 254 (184 up-regulated, 70 down-regulated) and 2,364 (1,247 up-regulated, 1,117 down-regulated) DEGs were detected in LT and Pdi, respectively. As it can be seen in Table 3, the largest proportion of genes in LT are annotated in cellular component organization or biogenesis, and the maximal percentage of genes in Pdi are annotated in catalytic activity. The majority of DEGs were up-regulated in both treatments (Figures 2A,B), with nearly 72.44 and 52.75% percentage of all the DEGs, respectively. In addition, the largest number of DEGs in LT was found at 48 h, and Pdi was at 24 h. Genes differentially expressed at two or more time points were defined as CEGs (co-differentially expressed genes), a total of six CEGs (three up-regulated genes and three down-regulated genes) appeared in both LT and Pdi at the four time points (Figure 2C), and five CEGs were annotated as encoding hypothetical protein, with only one gene (Kap003732) annotated as encoding geranylgeranyl transferase type-2 subunit beta. Here, we found the three up-regulated genes and the three down-regulated genes might be involved in biotic and abiotic stresses responses of K. apiculata.

Even though the counts of DEGs generated from qRT-PCR were not the same as FPKMs generated from RNA-seq, they nearly all showed the same trends (Figure 3). The high correlation between qRT-PCR and RNA-seq data were a validation of the generated RNA-seq results.

A total of 451 GO terms were obtained in the GO functional category analyses, the number of GO terms of biological process, cellular component and molecular function were 159, 161, and 131, respectively (Figure 4A). For biological process annotation, cellular process (GO: 0009987) and metabolic process (GO: 0008152) were the most plentiful terms, and terms annotated as response to stimulus (GO: 0050896) were also abundant. Meanwhile, binding (GO: 0005488), catalytic activity (GO: 0003824), and molecular function regulator (GO: 0098772) were largely associated terms in the annotation of molecular function. Moreover, cell part (GO: 0044464), protein-containing complex (GO: 0032991), and membrane-enclosed lumen (GO: 0031974) were the most involved terms for cellular component.

For KEGG enrichment analyses, ribosome biogenesis in eukaryotes (KEGG: sce03008), ABC transporters (KEGG: sce02010), and arginine biosynthesis (KEGG: sce00220) were the significantly enriched pathways, they all showed a higher enrichment factor (Figure 4B). In addition, phenylalanine metabolism (KEGG: sce00360) pathway was also enriched in the analyses.

There were 7 CEGs in the four time points in LT (Figure 2A), and the normalized expression counts of the 7 CEGs were visualized in a heatmap, and the expression patterns were clustered in STEM (Figures 4C,D), the 7 CEGs were grouped in 3 clusters (7, 4, 17), gene Kap003732 remained up-regulated throughout the entire treatments. Interestingly, we also observed one CEG (Kap003067) with annotation as heat shock protein 78 (HSP78) was downregulated at 5, 12, and 48 h, but it was not detected differentially expressed at 24 h.

Because there were many DEGs at different time points, the CEGs from Pdi were utilized for GO category to investigate the possible functions analyses. Results from GO category analyses were on a same occasion to LT, the CEGs were primarily annotated in biological process and cell components (Figure 5A). To get a better understanding of the CEGs, the entire 1,167 CEGs at the four time points in Pdi were employed for further analyses. The normalized expression counts of the 1,167 CEGs were visualized in a heatmap and expression patterns were clustered in STEM (Figures 5B,D), and the 1,167 CEGs were clustered into 20 different groups, group 17 (79 CEGs) was observed up-regulated throughout the entire treatments (Figure 5D). We also performed KEGG enrichment analyses to detect the potential biological pathways, and the results showed that many DEGs were annotated to metabolic pathways (KEGG: sce01100) and biosynthesis of secondary metabolites (sce01110) (Figure 5C). Moreover, lysine biosynthesis (KEGG: sce00300) pathway was enriched with the highest factor, and steroid biosynthesis (KEGG: sce00100) followed. CEGs involved in lysine biosynthesis (KEGG: sce00300), steroid biosynthesis (KEGG: sce00100), phenylalanine metabolism (KEGG: sce00360), and riboflavin metabolism (KEGG: sce00740) KEGG pathways were selected and the FPKM were normalized and visualized in a heatmap (Figure 6).

We also performed GO category and KEGG enrichment analyses of each individual time point. DEGs were associated with catalytic activity (GO: 0003824), hydrolase activity (GO: 0016787), and phosphatase activity (GO: 0016791) after co-incubation for 5, 12, and 24 h (Supplementary Figures 1–3). In addition, DEGs were involved in co-translational protein targeting to membrane (GO: 0006613), sodium ion transport (GO: 0006814), and others after co-incubation for 48 h (Supplementary Figure 4).

In the present study, two different kinds of treatments (or stresses) were carried out on K. apiculata yeast, and the expression levels were obtained by RNA-seq. A comparison was conducted using WEGO to extend a comprehensive understanding of these treatments, and the detailed GO terms were shown in Table 2. The number of DEGs annotated to catalytic activity (GO: 0003824) of LT and Pdi were 64 (25.4%) and 717 (32.1%), respectively. However, a higher percentage of DEGs annotated to catalytic activity (acting on RNA) (GO: 0140097) of LT than Pdi, with proportion of 6.6 and 3.9%, respectively. DEGs annotated to oxidoreductase activity (GO: 0016491) of LT and Pdi were 6 (2.4%) and 121 (5.4%), respectively. The percentage of DEGs annotated to the cellular component term of LT were mainly reduced compared to Pdi, excluding the endomembrane system.