Strain VD991 of V. dahliae was grown on potato dextrose agar (PDA) until sufficient spores were produced for sample collection. DNA was extracted from the spores using a cetyl-trimethyl ammonium bromide (CTAB) method (Biel and Parrish, 1986; Wang et al., 2022). The DNA obtained was then quality tested (Nanodrop and 0.35% agarose gel electrophoresis) and quantified (Qubit). Large fragments of genomic DNA were recovered via the BluePippin automatic nucleic acid recovery system, and DNA libraries were prepared using the SQK-LSK109 Ligation Sequencing Kit. Sequencing was performed using the PromethION Flow Cell (R9 Version) (Oxford Nanopore Technologies) (Deamer et al., 2016; Lu et al., 2016).

The assembly of the V. dahliae VD991 was carried out using the following process: longer reads and high-quality reads were extracted using Filtlong v0.2.0²; the filtered Nanopore reads were assembled using NECAT v0.01 (Chen et al., 2021); quality control of Illumina short reads was conducted using Trimmomatic v0.30 (Bolger et al., 2014); the reads obtained were used for polishing of the sequences assembled from Nanopore reads; six rounds of assembly polishing were carried out on Illumina reads using Pilon v1.23 (Walker et al., 2014) to correct base-calling and insertion/deletion errors. Hi-C fragment libraries were constructed and Illumina HiSeq sequencing was performed. Clean reads were mapped to the V. dahliae genome using BWA v0.7.9 (Li and Durbin, 2009). Paired-end reads were mapped to the genome separately and filtered, followed by the collection of unique, mapped paired-end reads using HiC-Pro v2.10 (Servant et al., 2015). The order and direction of scaffolds/contigs were clustered into super scaffolds using LACHESIS, based on the relationships among valid reads. Finally, the data were assembled onto eight chromosomes. By aligning next-generation sequencing data for V. dahliae VD991 against the fully assembled genome, genome quality was assessed based on the percentage and coverage of mapped reads. In addition, the BUSCO (Benchmarking Universal Single-Copy Orthologs) method (based on 290 conserved core genes for fungi) was used to further assess assembly completeness and quality (Simao et al., 2015).

As repetitive sequences tend to be poorly conserved among species, a species-specific repeat database was developed for V. dahliae VD991 using de novo and structural prediction as implemented in LTR_FINDER v1.05 (Xu and Wang, 2007), MITE-Hunter (Han and Wessler, 2010), PILER-DF v2.4 (Smith et al., 2007), and RepeatScout v1.0.5 (Price et al., 2005); PASTEClassifier v2.0 (Wicker et al., 2007) was then used to classify the repeated elements in the database. The newly constructed database was merged with the Repbase database (Bao et al., 2015), and repeated sequences were predicted for the V. dahliae VD991 genome using RepeatMasker v4.0.6 (Tarailo-Graovac and Chen, 2009).

Prediction of gene structure was performed using de novo prediction, homologous protein-based prediction, and transcriptome-based prediction. De novo prediction was conducted using Augustus v2.4 (Stanke and Waack, 2003), GeneID v1.4 (Blanco et al., 2007), Genscan (Burge and Karlin, 1997), GlimmerHMM v3.0.4 (Majoros et al., 2004), and SNAP v2006-07-28 (Korf, 2004), while homologous protein-based prediction was implemented in GeMoMa v1.3.1 (Keilwagen et al., 2016). The transcriptomic data were assembled using Hisat2 v2.0.4 (Pertea et al., 2016) and Stringtie v1.2.3 (Pertea et al., 2016), and unigene sequences predicted using PASA v2.0.2 (Campbell et al., 2006) and TransDecoder v2.0 (Tian et al., 2018). Finally, EVM v1.1.1 (Haas et al., 2008) was used to integrate the results from all three prediction methods.

Non-coding RNAs are RNAs that do not encode proteins, including RNAs with known functions such as microRNAs, rRNAs, and tRNAs. Using known structural characteristics of non-coding RNAs, tRNAs were predicted for the V. dahliae genome using tRNAscan-SE (Lowe and Eddy, 1997); rRNAs and other ncRNAs (i.e., not rRNAs and tRNAs) were predicted using Infernal v1.1 (Nawrocki and Eddy, 2013) and the Rfam database (Nawrocki et al., 2015).

Pseudogenes are similar in sequence to functional genes, but have become nonfunctional due to mutations such as insertions and deletions. To identify homologous gene sequences in the V. dahliae VD991 genome, predicted protein sequences were compared with protein sequences from the Swiss-Prot database using GenBlastA v1.0.4 (She et al., 2009). To identify pseudogenes, GeneWise (Birney et al., 2004) was used to find premature stop codons and frameshift mutations in known gene sequences.

The predicted gene sequences were aligned to functional databases [e.g., KEGG (Kanehisa et al., 2004), KOG (Tatusov et al., 2000), Nr (Deng et al., 2006), Swiss-Prot (Boeckmann et al., 2003) and TrEMBL (Boeckmann et al., 2003)] using BLAST v2.2.26 (Altschul et al., 1997) to annotate functional genes. Based on the Nr database annotation results, gene ontology (GO) annotation was performed using Blast2GO v5.2.5 (Conesa et al., 2005) and the GO database (Ashburner et al., 2000), while Pfam annotation was performed using hmmer v3.3.2 (Eddy, 1998) and the Pfam database (Finn et al., 2016). In addition, GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway enrichment analyses were also performed.

Protein sequences for predicted genes were annotated via alignment to functional databases such as the Pathogen-Host Interaction Factor Database (PHI) (Winnenburg et al., 2006) and the Transporter Taxonomy Database (TCDB) (Saier et al., 2006). In addition, functional annotation of carbohydrate-related enzymes based on the Carbohydrate-Associated Enzyme Database (CAZy) (Cantarel et al., 2009) was performed using hmmer v3.3.2 (Eddy, 1998).

Signal peptides are short peptide chains (usually 5–30 amino acids in length) that guide the transfer of newly synthesized proteins to the secretory pathway. The protein sequences of all predicted genes were analyzed using SignalP v4.0 (Petersen et al., 2011) and tmhmm v2.0 (Krogh et al., 2001) to identify proteins containing signal peptides and/or transmembrane helices (i.e., transmembrane proteins), respectively. To verify the accuracy of the annotated genes obtained by bioinformatics prediction, two highly expressed genes that are both secreted and effector proteins were selected for sequencing. We designed primers (Supplementary Table 2) based on the gene sequences, extracted RNA from V. dahliae and reverse transcribed it into cDNA. cDNA was then used as a template for gene amplification and sanger sequencing.

Three varieties of cotton with varying degrees of disease resistance were selected for transcriptomic analysis: the susceptible variety Jimian11 (JM11), the resistant variety Zhongzhimian2 (ZZM2), and Zhongmiansuo24 (ZM24). Seedlings from each strain were grown in the greenhouse for 15 days and transferred into V. dahliae broth upon emergence of the first true leaf (Li et al., 2021). Seedling roots were sampled at 0, 6, 12, and 24 h after transfer into V. dahliae broth (The time point of 0 h represented that the roots had not yet been inoculated with V. dahliae, and we sampled and sequenced the roots and V. dahliae separately). Root samples from two seedlings were combined to form a single sample, and there were three biological replicates for each strain. Root samples were sent to Beijing Biomarker Company for transcriptome sequencing.

The previously published TM-1 genome (Yang et al., 2019) and the newly assembled V. dahliae VD991 genome (obtained in this study) were used as reference genome. Transcriptome sequencing data were mapped to the reference genome, and count matrices and FPKM (Fragments Per Kilobase Million) values were obtained using previously described methods (Ding et al., 2019). The values from roots and V. dahliae separately at 0 h were combined as a control group. Differentially expressed genes (DEGs) were identified as those with a P-value < 0.05 and | log2-FoldChange| > 1.2. Gene ontology (Ashburner et al., 2000) and KEGG (Kanehisa et al., 2004) enrichment analyses were performed using the R-package “clusterProfiler”. In addition, a weighted gene co-expression network analysis (WGCNA) was performed using the R package “WGCNA” (Langfelder and Horvath, 2008; Duan et al., 2022), with the FPKM values as input, as described previously (Schurack et al., 2021).

Eighty-seven V. dahliae resequencing datasets were downloaded from the NCBI SRA database³; the V. dahliae VD991 genome, newly assembled in this study, was then used as the reference genome for resequencing analysis (Yang et al., 2019). After calling SNPs, samples were grouped based on the results of a principal components analysis (PCA), population structure analysis and phylogenetic analysis. The population fixation index (F_ST) was calculated (based on the grouping results) using vcftools v0.1.13 (Danecek et al., 2011).

A total of 5.09 Gb of raw reads were obtained from Nanopore sequencing of V. dahliae VD991. After removing adapters, short fragments, and low-quality data, a total of 4.92 Gb (∼137.4 ×) of clean reads were obtained for use in whole genome assembly. The final assembled genome for V. dahliae VD991 consisted of nine scaffolds and a scaffold N50 length of 4,119,679 bp, with the longest scaffold having a length of 7,830,508 bp and a GC content of 53.41%. Roughly 35.77 Mb of sequence data were anchored onto eight pseudochromosomes, with 99.92% of the sequences oriented (Table 1 and Figure 2A).

Assembly completeness was assessed by alignment with Illumina reads and BUSCO analysis. More than 99.17% of Illumina reads mapped properly to the new assembly. Furthermore, 98.62% of 290 core conserved genes (from the fungi_odb9 database) were classified as complete in the BUSCO analysis.

Roughly a tenth (9.33%) of the assembled genome was classified as repetitive sequences (Supplementary Table 3). Class I transposable elements (TEs) constituted the predominant repeat type, accounting for 8.77% of the total genome length. A total of 10,455 genes were predicted (Supplementary Table 4) by combined de novo, homologous protein-based, and transcriptome-based prediction; of these, 10,441 (99.86%) genes were supported by both homologous protein-based and transcriptome-based predictions, suggesting these are well-supported genes (Supplementary Figure 1). The average gene length was 2,142 bp, with an average of 2.94 exons, 1.94 introns, and 2.87 CDS per gene. Prediction results for non-coding RNA identified 125 rRNAs, 247 tRNAs, and 36 additional unclassified non-coding RNAs. The secreted and effector proteins are considered important pathogenic factors of V. dahliae. In this study, 854 secreted proteins (Supplementary Table 5) and 128 effector proteins (Supplementary Table 6) were predicted. To validate the accuracy of the gene prediction, cDNA amplification and sequencing were performed on the predicted effector and secretory protein genes, and the amplified gene sequences were found to be consistent with the predicted sequences (Supplementary Figures 2, 3), confirming the accuracy of our gene prediction results.

Gene expression levels were significantly lower in the 10-kb region surrounding each TE vs. other regions (Figure 2B). To assess how TE type affected gene expression, genes located in the vicinity of TEs were divided into those near Class I TEs only and those near Class I and Class II TEs. Genes affected only by Class I TEs showed significantly lower expression as compared to genes affected by both Class I and Class II TEs (Figure 2C). Thus, Class I TEs might suppress gene expression, while Class II TEs counteract this effect (to some extent). Transposable elements are an important source of mutations and genetic polymorphism. Many TE families are still actively transposable, and this process is highly mutagenic. In animals, plants, and microorganisms, many mutations (and the resulting phenotypic variation) are caused by transposition of these elements (Bourque et al., 2018). By calculating the density of TEs near indels (Figure 2A), more TEs were found in the vicinity of indels vs. other regions, and the density of TEs increased as the distance to the indel decreased (Figure 2D).

The genomes of 25 closely related species to V. dahliae, including 18 Colletotrichum species, two Plectosphaerella (ascomycetes), one Sodiomyces and Verticillium fungicola (four genomes), were used to construct a phylogenetic tree. In the tree, V. dahliae diverges early within this evolutionary lineage. Computational analysis of gene family evolution (CAFE) was used to estimate the number of gene families that have experienced historical expansion or contraction; 53 gene families were found to have expanded and 285 gene families to have contracted in V. dahliae (Figure 3A). A collinearity analysis revealed that the genomes of Plectosphaerella cucumerina and Verticillium alfalfae partially overlapped with the genome of V. dahliae, suggesting that a whole genome duplication (WGD) event occurred in V. dahliae (Figure 3B). Searching for further genome duplication events, both V. dahliae and Verticillium alfalfae showed peaks at 4DTV = 0.05 (Figure 3C), and this finding was further supported by the Ka/Ks analysis (Figure 3D).

Using 87 V. dahliae genomes downloaded from the NCBI database, including both deciduous and non-deciduous types, pathogenicity-related genes were explored for V. dahliae: 302,949 single nucleotide polymorphisms (SNP) were identified. Phylogenetic and structural analyses based on the SNP data divided the 87 accessions into two subgroups (Figure 4A), and this division was further supported by the PCA (Figure 4B). Linkage disequilibrium (LD) analysis was used to quantify the genetic diversity within populations. Linkage disequilibrium decayed more slowly in the high-toxicity population (G2) than in the low-toxicity population (G1), indicating less genetic diversity in the high-toxicity population (Figure 4C).

Most SNPs were located in intergenic regions, suggesting they do not directly affect gene structure (Figure 4D). V. dahliae VD991 as a reference genome is a highly virulent strain. SNP density was measured for both subgroups (see above), and a lower SNP density was found in the high-toxicity population (Supplementary Figure 4), further reinforcing the accuracy of the grouping. To identify pathogenicity-related genes in V. dahliae, F_ST values were calculated between the high- and low-toxicity populations; selecting sites with the top 5% F_ST values, a significant region was located on chromosome 6 (Supplementary Figure 5). Annotation of this region revealed five potentially pathogenic genes, including a glucose/galactose transporter gene (Vd06G0684, Vd06G0688, Vd06G0691, Vd06G0693, and Vd06G0694) (Figure 4E and Supplementary Table 7). In addition, a KEGG analysis of all genes located within this region found significant enrichment of a number of pathways associated with disease resistance, such as yeast cell cycling and glycosaminoglycan degradation (Figure 4F).

To investigate the interactions between V. dahliae and cotton, three cotton varieties (a susceptible variety JM11, resistant variety ZZM2, and ZM24) were inoculated with V. dahliae VD991 and transcriptome sequencing was performed. A total of 153.46 Gb of raw sequencing data were obtained containing 59,730 genes with large changes in expression. A WGCNA was used to investigate gene expression (in the genes with large changes in expression) at different times after V. dahliae inoculation, resulting in twelve genetic modules (Figures 5A, B). Effector proteins are known to play an important role in V. dahliae infection (Stergiopoulos and Wit, 2009; Feng et al., 2018; Wang et al., 2020); therefore, the cyan module containing the most effector protein genes was selected for subsequent analysis. By examining genes in V. dahliae important for interactions with cotton, pairs of interacting genes were obtained with weights greater than 0.25. A hub gene mining analysis was then performed using the MCODE package in Cytoscape; this resulted in a network containing 19 hub genes that may interact with cotton genes (Figure 5C). GO and KEGG enrichment analyses were performed of all the cotton genes interacting with V. dahliae. A total of 216 significantly enriched GO entries were obtained using P < 0.01. The top 10 entries were selected (with the smallest P-values, as shown in Figure 5D). In addition, the top 10 KEGG pathways were selected with P-values < 0.01 and gene counts (Figure 5D). In the GO and KEGG analyses, several terms and pathways were associated with resistance to V. dahliae, including the response to oxidative stress (GO:0006979) and phenylalanine metabolism (ko00360).

The darkslateblue module was positively correlated with V. dahliae post-inoculation time points (0, 6, 12, and 24 h); therefore, 5,000 gene pairs (with the highest weights) were selected from this module for analysis, resulting in a network containing 176 hub genes (Supplementary Figure 6). Ubiquitination plays an important role in plant resistance to pathogen invasion (Gao et al., 2022; Li et al., 2022). Here, 10 cotton hub genes (Gh_A08G270200, Gh_A09G093700, Gh_A09G208700, Gh_A09G238700, Gh_A10G060400, Gh_A13G194800, Gh_A13G255500, Gh_D03G157000, Gh_D03G170700, and Gh_D13G260100) were associated with ubiquitination. Two cotton hub genes (Gh_A12G006900 and Gh_D12G006300) were related to autophagy. In addition, a number of cotton genes associated with disease resistance were identified, such as genes involved in the jasmonic acid pathway (Gh_A06G223900) and catalase hydrogen peroxide (Gh_D06G205800).

For each post-inoculation timepoint, DEGs were identified between JM11 and ZZM2; DEGs were also identified between adjacent time points for JM11 and ZZM2 (individually). The greatest number of DEGs between JM11 and ZZM2 occurred at 6 h post-inoculation, with more up-regulated DEGs than down-regulated DEGs. By 12 h and 24 h, the number of DEGs had declined, and there were fewer up-regulated DEGs vs. down-regulated DEGs (Figure 6A). In V. dahliae, the number of DEGs initially increased to 6 h post-inoculation, then decreased significantly by 12 h before again increasing to a maximum at 24 h. Most of the DEGs in V. dahliae were down-regulated at 6 h, with almost all of the DEGs being down-regulated by 12 and 24 h (Figure 6A and Supplementary Figure 7). The common genes at these three time points and 0 h were removed (Figure 6B), leaving the remaining DEGs associated with disease resistance. A KEGG enrichment analysis of the common DEGs identified the following enriched pathways: ubiquinone and other terpenoid quinone biosynthesis, phenylpropanoid biosynthesis, and phenylalanine metabolism (Figure 6C). One of the pathways associated with disease resistance (phenylalanine metabolism) was illustrated to show the expression of DEGs in this pathway. As shown in Figure 6D, six genes were differentially expressed between JM11 and ZZM2, suggesting they may underlie the variation in disease resistance among the cotton varieties.