Comparative Genomics Reveals Evolutionary Traits, Mating Strategies, and Pathogenicity-Related Genes Variation of Botryosphaeriaceae.

Botryosphaeriaceae, as a major family of the largest class of kingdom fungi Dothideomycetes, encompasses phytopathogens, saprobes, and endophytes. Many members of this family are opportunistic phytopathogens with a wide host range and worldwide geographical distribution, and can infect many economically important plants, including food crops and raw material plants for biofuel production. To date, however, little is known about the family evolutionary characterization, mating strategies, and pathogenicity-related genes variation from a comparative genome perspective. Here, we conducted a large-scale whole-genome comparison of 271 Dothideomycetes, including 19 species in Botryosphaeriaceae. The comparative genome analysis provided a clear classification of Botryosphaeriaceae in Dothideomycetes and indicated that the evolution of lifestyle within Dothideomycetes underwent four major transitions from non-phytopathogenic to phytopathogenic. Mating strategies analysis demonstrated that at least 3 transitions were found within Botryosphaeriaceae from heterothallism to homothallism. Additionally, pathogenicity-related genes contents in different genera varied greatly, indicative of genus-lineage expansion within Botryosphaeriaceae. These findings shed new light on evolutionary traits, mating strategies and pathogenicity-related genes variation of Botryosphaeriaceae.

Introduction

Dothideomycetes represents the largest and most important class of ascomycete fungi, including 23 orders, 110 families, 1,261 genera, and 19,000 species (Wijayawardene et al., 2017). The members of Dothideomycetes comprise both phytopathogenic (Ohm et al., 2012) and non-phytopathogenic fungi with diverse lifestyles (Ruibal et al., 2009) as well as many mycorrhizal fungi (Peter et al., 2016). Among these 110 families, Botryosphaeriaceae is an important and distinctive family. This is because it includes saprobes, endophytes, and phytopathogens, and it is one of the most widely geographically distributed groups of opportunistic plant pathogens. The host range of this family is very wide, and many economically important plants worldwide can be infected by them (Slippers and Wingfield, 2007). These pathogenic fungi can infect plants through wounds or natural openings, such as lenticels and stomata. Once they enter host tissues, they may survive as endophytes to stay at a biotrophic stage for a long time and turn into the destructive necrotrophic stage when the host is stressed (Yan et al., 2013; Morales-Cruz et al., 2015). The members of Botryosphaeriaceae can infect many woody plants and cause serious disease symptoms, such as dieback, branch canker, leaf spots, and fruit and seed rot (Marsberg et al., 2017). But the interaction of some of Botryosphaeriaceae fungi like Botryosphaeria dothidea (B. dothidea) with host plants includes a latent or endophytic phase, which makes the fungal infection easily be neglected (Bostock et al., 2014; Marsberg et al., 2017; Wang et al., 2018). Therefore, it is of substantial biological significance to explore the evolutionary characteristics of Botryosphaeriaceae at the level of Dothideomycetes.

The exceptional feature of Botryosphaeriaceae fungi is that it is difficult to observe their sexual structure under both natural and experimental conditions, but this does not mean that sexual reproduction does not occur for these fungi (Phillips et al., 2013). With more in-depth research being conducted, the mating strategies (homothallism or heterothallism) of increasing members of the Botryosphaeriaceae family have been reported (Bihon et al., 2012a,b, 2014; Billiard et al., 2012; Lopes et al., 2017, 2018), revealing some unresolved questions regarding their sexual reproduction. For example, how conservative are the nucleic acid and protein sequences of the mating type determination genes; how conservative are the genes and their arrangement at the mating type determination loci; how has the mating type evolved, and what is the origin type. In addition, the host range of different Botryosphaeriaceae fungi varies greatly. B. dothidea is a common pathogen with a wide range of hosts. Generally, the infection becomes symptomatic when the host is subjected to drought, physical damage, waterlogging, or freezing stress (Bostock et al., 2014; Marsberg et al., 2017; Wang et al., 2018). The symptoms primarily include canker on young seedlings, branches, and stems; the necrosis of branches; and fruit decay, which may lead to the death of the host in extreme cases (Kim et al., 2005; Tang et al., 2012). Botryosphaeria kuwatsukai (B. kuwatsukai) can also cause symptoms, such as fruit softening and decay, and severe canker of branches and stems (Wang et al., 2021). However, B. kuwatsukai has a relatively narrow host range and primarily infects apple and pear trees (Xu et al., 2015). Therefore, a systematic study of the mating strategies and differences of pathogenicity-related genes in Botryosphaeriaceae fungi will lead to a better understanding of their molecular evolutionary history and pathogenic characteristics.

Currently, genomics technology has been widely used to study many pathogenic fungi of plants, and greatly promoted the understanding of their evolution and pathogenic mechanisms (Islam et al., 2012; Blanco-Ulate et al., 2013; van der Nest et al., 2014; Nagel et al., 2018; Félix et al., 2019; Landi et al., 2020; Meile et al., 2020). For example, gene family expansion associated with virulence factors in wood-colonizing pathogenic fungi in the Botryosphaeriaceae was revealed via phylogenomic comparisons (Garcia et al., 2021); the comparative genome analyses of latent plant pathogens in the Botryosphaeriaceae were conducted to define their genomes (Nagel et al., 2021). In addition, a large-scale comparative genomic analysis can better reveal the physiological characteristics and evolutionary history of fungi (Guttman et al., 2014; Haridas et al., 2020; Miyauchi et al., 2020). Therefore, in this study, 167 Dothideomycetes fungi, including 19 Botryosphaeriaceae species, were fully sequenced, and a comparative genomics approach was used to comprehensively analyze the molecular evolution characteristics, mating strategies and pathogenicity within Botryosphaeriaceae fungi. The sequence differences of related genes benefit our understanding of the evolutionary history of Botryosphaeriaceae fungi and provide useful information for the prevention and control of the diseases caused by these fungi.

Materials and Methods

Fungal Strains

A total of 167 Dothideomycetes fungal strains, including 160 Botryosphaeriaceae fungi (19 species) (Supplementary Table 1), from the Fungal Strain Library of Shandong Agricultural University, Tai’an, China, were sequenced in this study. The sequencing data are found at NCBI (PRJNA777748).

Sequencing and Genome Assembly

The CTAB (hexadecyltrimethylammonium bromide) method (Murray and Thompson, 1980) was used to extract high-quality DNA. A 500 bp DNA fragment library was constructed, and PE150 sequencing was performed using Illumina HiSeq 4,000 (San Diego, CA, United States). The raw data obtained from sequencing were inputted into Trimmomatic v0.39 (Bolger et al., 2014) for quality control, and reads with an average quality of less than 30 were filtered. Jellyfish v2.3.0 (Marçais and Kingsford, 2011) was used to calculate k-mer distribution and GenomeScope v2.0 (Ranallo-Benavidez et al., 2020) was then used to assess genome size of each fungus. The assembly was carried out using SPAdes v3.13.1 (Bankevich et al., 2012).

Gene Prediction and Genome Annotation

RepeatModeler v2.0.1[1] was used to construct custome repeat libraries for each assembly, and RepeatMasker v4.1.1[2] was used to determine repeat contents with the custom repeat libraries. For ab initio gene prediction, GeneMark-ES v4.48_3.60 (Lomsadze et al., 2014) and AUGUSTUS v3.2.1 (Stanke and Morgenstern, 2005) were used with default parameters. First, GeneMark-ES was used to predict gene models, then the models were used to train AUGUSTUS. Exonerate v2.2.0 (Slater and Birney, 2005) was used for homology comparison prediction. Finally, MAKER v3.01.03 (Campbell et al., 2014) was used to predict protein-coding genes by combining the gene models from GeneMark-ES, AUGUSTUS, and Exonerate. BUSCO v4.1.4 was used to evaluate the completeness of genomes and genome annotations based on pezizomycotina_odb9 (3156 core ortholog genes) (Waterhouse et al., 2018). Functional annotations on the putative genes were performed using the following softwares: BLAST v2.8.1 + (Camacho et al., 2009) for the NCBI non-redundant protein (NR, 2020-06), SwissProt and FunSecKB2 databases (E-value threshold of 1E-05); HMMER v3.2.1 (Eddy, 2008) for the Pfam 32.0 and TransportDB 2.0 databases; KofamKOALA v1.2.0 (Aramaki et al., 2020) for KEGG annotation; and dbCAN2 (Zhang et al., 2018) for annotating carbohydrate active enzymes (CAZy).

Evolutionary Analysis of Botryosphaeriaceae Fungi

To construct an accurate evolutionary tree, eight outgroups (Supplementary Table 1) were selected and OrthoFinder v2.3.11 (Emms and Kelly, 2015) was used for gene family clustering. According to the clustering results, single copy orthogroups were extracted and aligned using MAFFT v7.427 (Katoh et al., 2002). These processed single copy orthogroups were then concatenated using a self-written Perl script and filtered using Gblock v0.91b (Castresana, 2000). PartitiosnFinder v2.1.1 (Lanfear et al., 2017) and RAxML v.8.2.2 (Stamatakis, 2014) were used to determine the optimal amino acid substitution model and build the Maximum Likelihood (ML) tree with 1,000 bootstrap replicates, respectively.

The known ecologies states (multistate and binary) and FUNGuild were used to determine the ecologies of each species (Nguyen et al., 2016). Mesquite v3.61 software was then used to infer ancestral ecological character states (Haridas et al., 2020).

Functional Enrichment Analysis of Botryosphaeriaceae

To find the differences in gene function annotation of Botryosphaeriaceae fungi with different lifestyles, scipy.stats package in Python was used to perform the two-tailed Fisher exact test for functional annotation with p-value threshold of 0.01. To reduce false positives, functional annotations with the number of genes less than 100 were filtered.

Mating Strategies of Botryosphaeriaceae

Mating Type Genes and Surrounding Genes

To determine the presence of mating type genes in the Botryosphaeriaceae fungal genome, the mating type genes of Diplodia sapinea (D. sapinea) and some neighboring genes (KF551229 and KF551228) (Bihon et al., 2012a) were used as templates to search homologous genes from the putative Botryosphaeriaceae fungal genes that has been generated in this study, using the partial alignment mode in BLASTx (Camacho et al., 2009). The mating type genes obtained were then inputted into the NCBI’s conserved domain database (Marchler-Bauer et al., 2015) to determine their functional domains.

Comparison of the Arrangement of Mating Type Loci

To compare the arrangement of mating type genes and their surrounding genes in the Botryosphaeriaceae fungal genome, a BLASTn alignment was conducted and then EasyFig version 2.2.2 (Sullivan et al., 2011) was used for collinearity analysis, with the E value threshold set to 1e–4.

Phylogenetic Comparison and Ancestral State Reconstruction

To study the mating strategies of Botryosphaeriaceae fungi, ML evolutionary trees of 24 Botryosphaeriaceae fungi were constructed with single-copy genes (using the same method described in section: Comparison of the Arrangement of Mating Type Loci). Briefly, OrthoFinder v2.3.11 (Emms and Kelly, 2015) was used to cluster gene families and extract single copy orthogroups. Each orthogroup was aligned and concatenated with MAFFT v7.427 (Katoh et al., 2002). Gblock v0.91b (Castresana, 2000) was used for filtering, and finally, RAxML v.8.2.2 (Stamatakis, 2014) was used to build an ML tree with the LG + I + G + F model. Mesquite v3.61[3] and the Mk1 likelihood model were used to reconstruct the evolutionary process of homothallism (hom) and heterothallism (het).

Changes in Genes Related to the Pathogenicity of Botryosphaeriaceae Fungi

The production of phytotoxic compounds in pathogenic fungi, such as secondary metabolite, secreted proteins, and carbohydrate-active enzymes, is one of the important infective weapons (Amselem et al., 2011). The genomes of Botryosphaeriaceae were searched for genes encoding the phytotoxic compounds and then the number difference of these genes was statistically analyzed as previously described (Wang et al., 2018).

Specifically, 23 Botryosphaeriaceae fungi (7 genera) and 10 other representative fungi were included to analyze changes in genes related to the pathogenicity in Botryosphaeriaceae. The 10 representative fungi contained 1 biotrophic fungus (Puccinia graminis), 2 necrotrophic fungi (Valsa mali and Pyrenophora triticirepentis), 2 saprophytic fungi (Neurospora crassa and Rhizopus oryzae), 3 hemibiotrophic fungi (Pyricularia oryzae, Colletotrichum higginsianum, and Zymoseptoria tritici), 1 symbiotic fungus (Laccaria bicolor), and 1 endophytic fungus (Peltaster fructicola).

Results

Genome Sequencing, Assemby, and Annotation

In this study, 167 genomes of Dothideomycetes, including 160 Botryosphaeriaceae (19 species), were sequenced, assembled, and annotated (Supplementary Table 1). All genomes were sequenced at average coverage 166 ± 48 x. The average assembled genome lengths of 167 Dothideomycetes ranged between 28.85 and 61.78 Mb, which were consistent with the sizes estimated by k-mer counting approach. The contig N50 values of assembled genomes varied from 41.81 to 779.02 kb with a mean of 241.77 kb. The repeat contents of assembled genomes varied from 0.95 to 10.98% with a mean of 5.18%. All assembled genomes have a high completeness with an average of 94.8 ± 4.3%, and the similar result was also found in genome annotations (average of 98.0 ± 2.9%).

Classification Based on Whole-Genome Data

To better understand the evolutionary characteristics of Botryosphaeriaceae, we used 167 newly sequenced in this study and 112 genomes available in public database to construct a whole-genome phylogenetic tree using 480 single-copy gene families. The phylogenetic tree was clearly divided into 2 clades, corresponding to 2 subclasses, namely Pleosporomycetidae and Dothideomycetidae (Supplementary Figure 1). A comprehensive phylogenetic analysis revealed that the members of Botryosphaeriaceae that belonged to Botryosphaeriales were added to Pleosporomycetidae (Figure 1 and Supplementary Table 2).

FIGURE 1

Whole-genome-based phylogenetic tree of 271 species from Dothideomycetes and 8 outgroups. All bootstrap values are 100% except for those shown. The orders of Dothideomycetes were well classified and were displayed by different colors. The two circles left of species names standed for lifestyle classification according to organism data and FunGuild, respectively.

Saprophytic Fungi Have a Larger Genome Size Than Phytopathogenic Fungi

The genome sizes of Dothideomycetes ranged from 17 (Piedraia hortae) to 177 Mbp (Cenococcum geophilum), and 7,896–34,881 protein-coding genes were detected (Supplementary Figures 2, 3 and Supplementary Table 1). Compared with non-pathogenic fungi, pathogenic fungi usually have a smaller genome (Supplementary Figure 2). The genome size of Botryosphaeriaceae is in the range of 28.85 (Aureobasidium pullulans) −61.78 Mbp (Dothiorella sarmentorum) (47.19 Mbp on average), and 11,505–16,851 proteion-coding genes were predicted (13,765 genes in average).

Saprophytic Fungi Are the Possible Evolutionary Ancestors of Botryosphaeriaceae Fungi

To infer the ecological characteristics of ancestors of Botryosphaeriaceae fungi, we analyzed their lifestyle evolution process at the class level. The results showed that the ancestral lifestyle was likely to be the saprophytic type (Figure 2), which is also supported by the maximum likelihood analysis (Supplementary Figure 4 and Supplementary Table 3). During the evolution of Dothideomycetes fungi, at least 6 transitions from non-phytopathogenic (NPP) to phytopathogenic (PP) were detected, including 4 major transitions, which are presented at the MRCA node of Mycosphaerellaceae (Dothideomycetidae) (Node 225: NPP = 0.1868, PP = 0.8132), Venturia (Pleosporomycetidae) (Node 189: NPP = 0.0039, PP = 0.9961) and Botryosphaeriales (Node 20: NPP = 0.0984, PP = 0.9016), along with the branching point of Setomelanomma-Bipolaris (Node 98: NPP = 0.0074, PP = 0.9926). In addition, the Botryosphaeriaceae fungi have undergone at least 3 transitions from saprophytic to pathogenic fungi, including Node 176 (NPP = 0.9859, PP = 0.0141), Node 180 (NPP = 0.9992, PP = 0.0008) and Node 246 (NPP = 0.9101, PP = 0.0899) (Supplementary Table 3).

FIGURE 2

Reconstruction of ancestral lifestyle character state of Dothideomycetes using Mesquite based on parsimony model as saprobe. Major lifestyle shifts were marked by six red star symbols.

Differences in the Gene Families of Plant Pathogenic and Saprophytic Fungi

To investigate the differences between plant pathogenic and saprophytic fungi from the perspective of gene family contractions and expansions, we used Fisher’s exact test to perform these analyses on gene families. Differences in 78 Pfam and 58 GO annotations were found between saprophytic and plant pathogenic fungi. In the saprophytic fungi, 42 Pfam and 35 GO terms showed more than 20% of expansion, and 31 Pfam and 14 GO terms displayed more than 20% of contraction. In plant pathogenic fungi, 32 Pfam and 17 GO terms showed greater than 20% of expansion, and 35 Pfam and 21 terms displayed more than 20% of contraction. Compared with saprophytic fungi, the plant pathogenic fungi contained more gene families that showed contraction. Compared with saprophytic fungi, 94 Pfam and 67 GO terms of Botryosphaeriaceae fungi showed greater than 20% of expansion, while 52 Pfam and 44 GO terms displayed contractions (> 20%) (Supplementary Table 4).

Mating Strategies of Botryosphaeriaceae Fungi

Mating Type Genes

To determine the mating type genes, MAT1-1 and MAT1-2 of 24 fungi (7 genera and 19 species) were analyzed in the Botryosphaeriaceae family. The results showed that the Botryosphaeria, Neofusicoccum, and Dothiorella genomes harbored both MAT1-1 and MAT1-2, indicating that they were homothallic. Diplodia, Macrophomina, and Neoscytalidium genomes harbored either MAT1-1 or MAT1-2, indicating that they were heterothallic fungi. Lasiodiplodia fungi had various mating strategies, including homothallism (L. gonubiensis) and heterothallism (Lasiodiplodia citricola, L. pseudotheobromae, and L. theobromae) (Supplementary Table 5). In addition, the protein domains of MAT1-1-1 and MAT1-2-1 were highly conservative. All the MAT1-1-1 proteins contain MATalpha domains, and all the MAT1-2-1 proteins contain MAT_HMG-box domains (Supplementary Table 6). Compared with the neighboring genes, the nucleic acid sequences of the mating type genes were poorly conserved (Supplementary Table 7), while the length of the coding sequence and the position and size of the intron were largely conservative.

Arrangement of Mating Loci

Arrangement analysis of mating loci showed that there were three types of arrangements at the mating type determining loci of Botryosphaeriaceae (Figure 3). In most Botryosphaeriaceae genomes (e.g., Diplodia sapinea), the MAT1 gene was primarily located between a collinear region that contains 4 protein-coding genes and 1 putative integral membrane (PIM) protein (Figure 3). The PIM contains pleckstrin homology domains and DUF2404. The four genes in the collinear region encode, in order, a DNA lyase (APN2), cytochrome c oxidase subunit VIa (CoxVIa), anaphase-promoting complex subunit 5 (APC5), and complex I intermediate-associated protein 30 (CIA30). According to the positional relation between the MAT1 and APN2 genes (the MAT1 gene was located at the proximal end of APN2), the second type of arrangement was observed in Botryosphaeria, L. gonubiensis (Lasiodiplodia), M. phaseolina (Macrophomina), and N. dimidiatum (Neoscytalidium). For example, the MAT1 gene of M. phaseolina and B. dothidea was located upstream of the APN2 gene, while in L. gonubiensis and N. dimidiatum, the MAT1 gene was located downstream of the APN2 gene, and the four genes in the collinear region were arranged in reverse order. However, the opposite arrangement was observed in the two strains of B. kuwatsukai-the MAT1 gene of the PG2 strain was located downstream of the APN2 gene, while LW030101 was located upstream. A similar type of reversed arrangement as the one found for B. kuwatsukai was also identified in two B. dothidea strains. The third type of arrangement was primarily observed for Neofusicoccum species in which the two subtypes of MAT1, MAT1-1, and MAT1-2, were not located in conjunction; they were either located both distantly from the collinear region or at different chromosomes (or scaffolds) (Figure 3).

FIGURE 3

Pairwise mating type and surrounding genes comparison between species of Botryosphaeriaceae. Genes (color coded arrows) were on genomic sequences (horizontal lines). Organization of genes were indicated by gray box. Abbreviations of genes: putative integral membrane protein containing DUF2404 domain (DUF2404), DNA lyase (APN2), cytochrome C oxidase subunit Via (Cox), anaphase-promoting complex subunit 5 (APC5) and Complex I intermediate-associated protein 30 (CIA30).

Reconstruction of the Ancestral State of the Mating Strategy

To understand the evolutionary characteristics of the mating type of Botryosphaeriaceae, we selected 24 representative Botryosphaeriaceae strains to reconstruct the evolutionary process of homothallism and heterothallism mating strategies in Botryosphaeriaceae (Supplementary Figure 5). Our analysis showed that het mating was likely to be the ancestral type (Node 2, het = 0.6198, hom = 0.3802) (Supplementary Figure 5 and Supplementary Table 8). During the course of evolution, Botryosphaeriaceae fungi have undergone at least 2 major transitions to their hom strategy. The first was located at the branching point of Lasiodiplodia theobromae and L. citricola (Node 14: het = 0.8291, hom = 0.1709), and the other transition was primarily observed at the branching point of Neofusicoccum and Dothiorella (Node 26: het = 0.4344, hom = 0.5656) (Supplementary Figure 5 and Supplementary Table 8).

Changes in Genes Related to the Pathogenicity of Botryosphaeriaceae

Secondary Metabolism

To find the key enzymes involved in the synthesis of secondary metabolites in Botryosphaeriaceae fungi, we used Pfam annotations to find 3 types of genes that encode these key enzymes, such as polyketide synthase (PKS), non-ribosomal peptide synthetase (NRPS), and dimethylallyl tryptophane synthase (DMATS). The PKS genes of Botryosphaeriaceae fungi differed significantly from each other (p = 2.0 × 10–16) (Figure 4A). Among all the species, the Macrophomina fungi contained the largest number of PKS genes (31 on average), while the Dothiorella fungi contained the least number (8 on average) (Supplementary Table 9). The genes related to secondary metabolites in fungi also included those that encode cytochrome P450 enzymes, regulatory factors, and transporters. Cytochromes P450 can catalyze the transformation of hydrophobic intermediates in the primary and secondary metabolic pathways and plays an important role in fungi. Compared with other fungi, Botryosphaeriaceae fungi contained the largest number of genes that encode cytochromes P450 (133-267). Among all the species, Neofusicoccum fungi contained the highest number of cytochromes P450 (267), followed by Botryosphaeria (n = 251) (Figure 4B). The ATP-binding cassette (ABC) or major facilitator superfamily (MFS) gene families played important roles in the transport of secondary metabolites. The numbers of ABC or MFS in Botryosphaeriaceae fungi were higher than those of other species (ABC: 268-353; MFS: 16-62). Among all the species, Botryosphaeria fungi contained the largest number of MFS (n = 62), while Neoscytalidium dimidiatum fungi contain the least amount (n = 16) (Figure 4B).

FIGURE 4

Comparing backbone and related genes of secondary metabolism among Botryosphaeriaceae fungi and other 10 fungal species. (A) Key gene family of secondary metabolism. In each column, Z-score was used to describe the trend of over-represented (+ 4–0) and down-represented (0 to −4) gene family. (B) Gene family number comparison of MFS_1, ABC_tran and CYPs among Botryosphaeriaceae fungi and other 10 fungal species.

Secreted Proteins

Pathogens can secrete a battery of proteins, which are deployed to the host-pathogen interface during infection, and these secreted proteins played important roles in fungal pathogenicity (O’Connell et al., 2012; Fernandes et al., 2014; Félix et al., 2016). In this study, we predicted the secreted proteins of Botryosphaeriaceae (Figure 5). The number of secreted proteins in Botryosphaeria (the B. dothidea fungi contain 1,034, and the B. kuwatsukai fungi contain 939) was not significantly different from Macrophomina (n = 1,118), and the B. kuwatsukai fungi seemed to secrete less proteins. The number of proteins secreted by the Neofusicoccum fungi varied greatly (922–1,028). Among all the species, the N. parvum fungi contained the largest number (n = 1,028), and N. cordaticola fungi contained the least (n = 922). Compared with other fungal genera, the Diplodia genus contained fewer secreted proteins (671–866) (Figure 5 and Supplementary Table 10).

FIGURE 5

Comparison of secreted proteins among Botryosphaeriaceae fungi and other 10 fungal species.

The transfer of secreted effector proteins to host plant cells was the key to pathogenesis of many plant pathogenic microorganisms. Secreted proteins less than 200 amino acids in length and rich in cysteine were considered as candidate secreted effectors (Wang et al., 2018). The Botryosphaeriaceae family contained two genera Macrophomina and Botryosphaeria that may release a large number of secreted proteins, and the number of small, secreted proteins of these two genera was also higher than that of other genera. For example, M. phaseolina (261), B. dothidea (253), and B. kuwatsukai (230) displayed significant expansions of gene families. Small secretory proteins rich in cysteine were the most abundant in Botryosphaeria fungi and the least in Diplodia fungi (79-113) (Figure 5 and Supplementary Table 10). In Botryosphaeria, both B. dothidea and B. kuwatsukai contained 154 small secretory proteins rich in cysteine, and the Pfam annotations of them included 19 and 21 known functional domains, respectively (Supplementary Table 11). By comparing the Pfam annotations of B. dothidea and B. kuwatsukai secretory proteins, we found that ribonuclease, the cerato-platanin family, and cysteine-rich secretory protein (CAP) family were only present in B. kuwatsukai, while the cell wall integrity and stress response component (WSC) domain, a putative carbohydrate binding domain was only present in B. dothidea (Supplementary Table 11).

Pathogenic and saprophytic fungi can secrete peptidases into their surroundings to degrade a variety of host-related proteases (Wang et al., 2018). This degradation mechanism has potential benefits in eliminating the activity of antifungal host proteins and providing nutrients. Compared with other fungi, Botryosphaeriaceae fungi contained a higher number of secretory peptidases. Among all the species, the M. phaseolina (Macrophomina) fungi contained the largest number of secretory peptidases (n = 298), followed by B. dothidea (n = 264) and B. kuwatsukai (n = 235), but both were similar to the semi-biotrophic phytopathogenic pathogens C. higginsianum (n = 233) and P. oryzae (n = 240). Diplodia fungi contain the lowest number of secreted proteins (n = 192; p = 0.8 × 10–3) (Figure 5 and Supplementary Table 10).

Carbohydrate Active Enzymes

The ability to degrade complex carbohydrates in plants is an important aspect of the lifestyle of phytopathogenic fungi (Wang et al., 2018). Compared with other species, the fungi in Botryosphaeriaceae family contained a higher number of carbohydrate active enzyme-related gene families, such as genes that encode glycoside hydrolases (GHs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activities (AAs), and carbohydrate binding modules (CBMs). Some representative species include N. parvum (753), B. dothidea (750), N. kwambonambiense (724), L. theobromae (719), N. ribis (716), N. umdonicola (716), and N. cordaticola (709). Within the Botryosphaeriaceae family, N. parvum has the largest number of carbohydrate active enzyme-related gene families (753), followed by B. dothidea (750), and B. kuwatsukai (675). Most Diplodia fungi contain fewer carbohydrate activity enzyme-related gene families, such as D. scrobiculata (522), D. corticola (548), D. mutila (577), and D. sapinea (556). However, for some species, such as D. seriata (662), there was a small increase in the number (Figure 6A and Supplementary Table 12).

FIGURE 6

Comparing CAZy among Botryosphaeriaceae fungi and other 10 fungal species. (A) Six CAZy classes: CBMs, Carbohydrate-Binding Modules; PLs, Polysaccharide Lyases; GTs, Glycosyl Transferases; CEs, Carbohydrate Esterases; AAs, Auxiliary Activities; GHs, Glycoside Hydrolases. (B) Distribution of CAZy related to cellulose, hemicellulose and lignin degradation. (C) Comparison of selected enzymes involved in PCW (plant cell wall) degradation among Botryosphaeriaceae fungi and other 10 fungal species. In heatmaps, Z-score were used to describe over-represented (+ 4–0) and down-represented (0 to −4) gene family.

The plant cell wall is a complex network structure composed of different polysaccharides, including cellulose, hemicellulose, pectin, and lignin (Terrett and Dupree, 2019). Proteins encoded by fungal carbohydrate active genes and related auxiliary genes can degrade plant cell walls into simple monomers that are absorbed as the carbon source to provide energy for fungi (Islam et al., 2012; Wang et al., 2018). Compared with other fungi, those of Botryosphaeriaceae contain expanded gene families involved in the degradation of lignin, cellulose, hemicellulose, and pectin. For example, the gene family involved in the degradation of lignin (107 genes in average) was significantly larger than that in other species (54 genes in average; p = 2.0 × 10–4, t-test). The number of gene families involved in the degradation of hemicellulose and pectin also expanded significantly, with 67 vs. 34 (p = 2.0 × 10–3) and 70 vs. 32 (p = 2.0 × 10–3), respectively. Among the Botryosphaeriaceae fungi, compared with other genera, it was found that Neofusicoccum fungi contain the largest number of gene families involved in the degradation of plant cell walls (287–329), followed by Botryosphaeria (219–314) and Macrophomina (291–306), while Diplodia fungi contain the least number of gene families (215–259) (Figures 6B,C and Supplementary Table 13).

Discussion

In this study, we conducted a large-scale whole-genome sequencing of 167 Dothideomycetes fungi, including 19 species, and comprehensively analyzed evolutionary traits, mating strategies and changes in pathogenic genes in Dothideomycetes fungi. Our results can advance our understanding of the evolutionary history of Botryosphaeriaceae fungi.

Our results also confirmed that the 167 Dothideomycetes fungi can be divided into two subclasses, Pleosporomycetidae and Dothideomycetidae, and Botryosphaeriaceae belongs to Pleosporomycetidae, indicating that it is accurate to classify fungi using phylogenies based on phylogenomics (Haridas et al., 2020). Our results inferred the ancestral lifestyle of Botryosphaeriaceae fungus as saprobe, and these fungi have undergone at least three transitions from saprophytic to phytopathogenic states. This result is consistent with Schoch et al. (2009). They inferred that the Dothideomycetes fungi have experienced multiple transitions from saprophytic pathogens to lichens to phytopathogens, along with multiple transitions from terrestrial to aquatic lifestyles (Schoch et al., 2009). In addition, a larger genome size and a higher number of protein-coding genes are usually associated with saprophytic fungi compared with phytopathogenic fungi. This result is consistent with the findings of Schuelke et al. (2017). They found that in the Geosmithia genus, compared with non-pathogenic fungi, the pathogenic fungi in this family have smaller genomes (Schuelke et al., 2017).

Previous studies have shown that phytopathogenic fungi have a smaller genome and contain fewer protein-coding genes compared with saprophytic fungi, which is related to the expansion and contraction of gene families (Ohm et al., 2012; Haridas et al., 2018). In this study, many gene families of phytopathogenic fungi showed contractions, primarily including genes which contain the Pfam domain tetratricopeptide repeat (TPR). TPR can interact with a variety of proteins, such as the anaphase-promoting complex, NADPH oxidases, and HSP90-binding proteins (Das et al., 1998; Kaneko et al., 2006; Kondo et al., 2010). The TPR protein is also part of the plant hormone signaling pathways (Schapire et al., 2006). The family enrichment results found that multiple members of the TRP family proteins (TPR 1, 3, 8, 10, 11, 12, 16, and 19) contracted in phytopathogenic fungi. This contraction may be due to a reduction of signal related TPR proteins in pathogenic fungi. This is because during the process of pathogenic fungal infection, fungal signal related TPR proteins will be affected by plant hormone signals (Haridas et al., 2020). Some other Pfam domains also contracted significantly in phytopathogenic fungi, including PF17111 (the fungal N-terminal domain of STAND proteins), PF05729 (the NACHT domain), and PF14479 (involved in prion inhibition and propagation). These domains are in heterokaryons and play an important role in incompatibility (Paoletti and Saupe, 2009; Greenwald et al., 2010; Daskalov et al., 2012). The contraction of Pfam domains associated with incompatibility in the heterokaryons suggests that other strategies to reduce the level of incompatibility in heterokaryons also exist to improve the adaptability of pathogenic fungi (Ishikawa et al., 2012). Simultaneously, compared with saprophytic fungi, many GO terms are contracted in phytopathogenic fungi, including protein kinases (GO: 0004672), transcription factors (GO: 0003700), zinc ion binding (GO: 0008270), and regulation of nitrogen utilization (GO: 0006808). Given the complexity of amino acid biosynthetic pathways and energy requirements, fungi rely on the absorption of plant amino acids to conserve their own energy. This could be one of the reasons for the contraction of GO terms in plant pathogens (Staats et al., 2013; Mur et al., 2017).

Botryosphaeriaceae fungi display both types of mating strategies, homothallism and heterothallism, which are determined by the single locus MAT1 (Mur et al., 2017). In this study, the mating type genes of 24 Botryosphaeriaceae fungal isolates (19 species) were analyzed, and we found that these genes were not highly conserved in terms of their nucleotide and amino acid sequences. However, for two subtype genes, two domains, MATalpha_HMGbox (MAT1-1-1) and MATA_HMG-box (MAT1-2-1), are more conservative. This is consistent with the fact that the mating type genes of fungi have almost no differences within a species, but they are highly divergent between species (Turgeon, 1998). In many ascomycete fungi, it is a conservative trend that the MATalpha and MATA_HMG-box regions contain introns (Arie et al., 2000; Paoletti et al., 2005; Duong et al., 2013; de miccolis et al., 2016). In D. sapinea and D. seriata, the introns in the MATA_HMG-box region of the MAT1-2-1 gene were lost, but the amino acid sequences that flank the lost intron site remain intact (Nagel et al., 2018). This phenomenon is consistent with the intron loss model derived from Poly-A primed mRNA (Sharpton et al., 2008). In this study, with the exception of Neofusicoccum, the arrangement of MAT1-1 and MAT1-2 genes is highly conserved in six other genera. They all are located between a collinear region that contains four protein-coding genes and one PIM protein, which is the same arrangement as found in most Phyllostictaceae fungi (Bihon et al., 2014; Wang et al., 2016). The arrangement of mating type genes of Neofusicoccum fungi is exceptional. The locations of the two subtypes of mating type genes are not in conjunction; they are either both located distantly from the collinear region or at different chromosomes (or scaffolds) (Lopes et al., 2017). Similarly, this discontinuous arrangement of mating type genes also exists in Aspergillus nidulans (Galagan et al., 2005), Curvularia cymbopogonis (Galagan et al., 2005), and Neosartorya fischeri (Rydholm et al., 2007).

The mating type gene MAT1 has two subtypes—MAT1-1 and MAT1-2. When two subtypes are both harbored in the genome, the mating strategy is hom, and if only one is present, it is het (Idnurm, 2011). The repetitive sequence-mediated deletion of one or more mating type genes can cause non-directional changes in mating types, such as the transition from self-fertility to self-sterility (Wilken et al., 2014; Xu et al., 2016; Yun et al., 2017). Such repetitive sequences were not observed in this study. Therefore, Botryosphaeriaceae fungi are unlikely to undergo a non-directional change regarding mating types (Nagel et al., 2018). In addition, to further understand the evolutionary characteristics of the fungal mating type in Botryosphaeriaceae, we reconstructed the ancestral state of the Botryosphaeriaceae fungal mating type and found that the het mating strategy is the ancestral type. Moreover, the fungi in this family experienced a number of transitions to the homothallism strategy, a shift that is common in ascomycete fungi (Inderbitzin et al., 2005; Yokoyama et al., 2006; Nygren et al., 2011; Gioti et al., 2012).

In this study, many gene families of the Botryosphaeriaceae fungi have shown significant expansions and contractions, and this change is conducive to the adaptation of fungi to the living environment (Alkan et al., 2013). These contracted gene families include genes that encode secondary metabolite synthases, secreted proteins, and carbohydrate active enzymes. In the family of genes that encode secondary metabolite synthases, NPRS and PKS gene clusters are responsible for the synthesis of toxic peptides and the production of naphthalenone pentaketides, respectively, in Botryosphaeriaceae fungi (D. seriata, L. theobromae and N. parvum) and other pathogenic fungi (A. fumigatus, Diaporthe ampelina, Phaeomoniella chlamydospora, and Togninia minima) (Andolfi et al., 2011; Paolinelli-Alfonso et al., 2016). Although the gene clusters of secondary metabolites are not regulated during infection, a large number of products of these gene clusters, which are significantly expanded in the genome, may be involved in the induction of disease symptoms and host adaptation (Yan et al., 2018). Cytochrome P450 is a superfamily of monooxygenases. In addition to participating in the post-synthesis modification of a variety of metabolites, it can also promote the adaptation of fungi to specific ecological niches by altering potentially harmful chemicals in the environment (Siewers et al., 2005). Here, we found the expansion of this gene family, which can explain the wide host range of Botryosphaeriaceae fungi. This is because these genes are also associated with some physiological characters; thus, their expansions are likely to promote pathogenic evolution (Nierman et al., 2005). In this study, expansions are obvious for the genes that encode ABC transporters in Botryosphaeriaceae fungi, indicating that these fungi have evolved stronger virulence and capacity against plant defense compounds (Han et al., 2001; Luini et al., 2010). Sugar synergistic transporters belong to the MFS family, and they play important roles in fungal spore formation, intercellular communication, and pathogenicity (Gaur et al., 2008). In the Botryosphaeriaceae fungi, we found that the sugar synergistic transporter gene family showed different degrees of expansion, which probably makes it more adaptive to different host environments (e.g., different pH values), and more resistant when interacting with various plant pathogens (Zhang et al., 2012; Yin et al., 2015).

Previous studies have showed that secreted proteins play important roles in the infection process of pathogenic fungi (Cobos et al., 2010; O’Connell et al., 2012; Fernandes et al., 2014; Félix et al., 2016, 2019). In this study, the secretory protein gene families of Botryosphaeriaceae fungi have expanded to different degrees, but there are large differences between different genera, which may be related to the different infection ranges of Botryosphaeriaceae fungi (Wang et al., 2018). In addition, the gene families of carbohydrate active enzymes in Botryosphaeriaceae fungi also showed different degrees of expansion. Among these families, the glycoside hydrolase family GH33 is composed of sialidase, which can hydrolyze the glycosidic bonds of terminal sialic acid residues in oligosaccharides. Sialidase can function as a pathogenic factor, facilitating the adaption to the host by evading host recognition or inhibiting host defense responses (Alviano et al., 2004). The cell wall of most dicotyledonous plants is composed of approximately 35% pectin. Pectin-degrading enzymes contribute to the degradation of the cell wall. This local degradation of the cell wall is necessary for fungi to enter the plant cytoplasm and replicate in the host tissue (ten Have et al., 1998). Similar to the highly pathogenic Colletotrichum higginsianum, Neofusicoccum, and Botryosphaeria both possess a larger number of genes encoding pectin-degrading enzymes, which may explain the difference in pathogenicity between Botryosphaeriaceae fungi (Wang et al., 2018). Many genes encoding cellulase (AA9 and GH12) and hemicellulase (GH31 and GH43) have been significantly expanded in Neofusicoccum and Botryosphaeria fungi. These expansions may explain the rapid infection and colonization of Botryosphaeriaceae fungi in woody plants (Yan et al., 2018).

Conclusion

In conclusion, we constructed a phylogenetic tree using whole-genome data and clarified the taxonomic position of Botryosphaeriaceae in Dothideomycetes. Heterothallism is the ancestral mating type of Botryosphaeriaceae fungi, and these fungi have undergone at least 3 transitions from heterothallism to homothallism. The host range of Botryosphaeriaceae infection is closely related to the changes in the number of pathogenic genes. Our results provide important insights into the evolutionary history, mating strategies and pathogenicity-related genes variation in Botryosphaeriaceae.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://zenodo.org/record/5184447#.YRUFTci0yAc; the sequencing data are found at NCBI (PRJNA777748).

Author Contributions

CY and HL developed the concept of this study and were main contributors to writing the manuscript. CY, YD, QL, and JZ performed all experiments, data analysis, and prepared figures. CY, HL, SC, XX, AL, and XZ contributed to the manuscript edit and review. All authors read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.