Bei B. Lü1†
Guo G. Wu1†
Yu Sun1†
Liang S. Zhang2
Xiao Wu1Wei Jiang1
Peng Li1Yan N. Huang1Jin B. Wang1
Yong C. Zhao3Hua Liu1Li L. Song1Qin Mo1
Ai H. Pan1Yan Yang4Xuan Q. Long5Wei D. Cui5Chao Zhang6Xu Wang7,8Xue M. Tang1*- 1Biotechnology Research Institute, Key Laboratory of Agricultural Genetics and Breeding, Shanghai Academy of Agricultural Sciences, Shanghai, China
- 2Institute of Agriculture and Biotechnology, Zhejiang University, Hangzhou, China
- 3Institute of Edible Fungi, Yunnan Academy of Agricultural Sciences, Yunnan, China
- 4Institute of Edible Fungi, Shanghai Academy of Agricultural Sciences, Shanghai, China
- 5Institute of Microbiology, Xinjiang Academy of Agricultural Sciences, Ürümqi, China
- 6Translational Medical Center for Stem Cell Therapy and Institute for Regenerative Medicine, Shanghai East Hospital, School of Life Sciences and Technology, Tongji University, Shanghai, China
- 7Department of Pathobiology, Auburn University, Auburn, AL, United States
- 8HudsonAlpha Institute for Biotechnology, Huntsville, AL, United States
The precious rare edible fungus Morchella conica is popular worldwide for its rich nutrition, savory flavor, and varieties of bioactive components. Due to its high commercial, nutritional, and medicinal value, it has always been a hot spot. However, the molecular mechanism and endophytic bacterial communities in M. conica were poorly understood. In this study, we sequenced, assembled, and analyzed the genome of M. conica SH. Transcriptome analysis reveals significant differences between the mycelia and fruiting body. As shown in this study, 1,329 and 2,796 genes were specifically expressed in the mycelia and fruiting body, respectively. The Gene Ontology (GO) enrichment showed that RNA polymerase II transcription activity-related genes were enriched in the mycelium-specific gene cluster, and nucleotide binding-related genes were enriched in the fruiting body-specific gene cluster. Further analysis of differentially expressed genes in different development stages resulted in finding two groups with distinct expression patterns. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment displays that glycan degradation and ABC transporters were enriched in the group 1 with low expressed level in the mycelia, while taurine and hypotaurine metabolismand tyrosine metabolism-related genes were significantly enriched in the group 2 with high expressed level in mycelia. Moreover, a dynamic shift of bacterial communities in the developing fruiting body was detected by 16S rRNA sequencing, and co-expression analysis suggested that bacterial communities might play an important role in regulating gene expression. Taken together, our study provided a better understanding of the molecular biology of M. conica SH and direction for future research on artificial cultivation.
Introduction
Medicinal mushrooms have become an attractive option for hygienic food or as a source of developing pharmaceuticals and nutraceuticals (Wasser, 2017). Medicinal properties for human health may be due to the various cellular components and secondary metabolites of the fruiting body (Martin et al., 2010; Chen et al., 2012; Lo et al., 2012). True morels (Morchella spp.) are ascomycetous fungi with high reputation for their edibility and appearance, which is similar to “a sponge on a stick” (Vieira et al., 2016). Taxonomically, Morchella is a well-defined genus in the class of Pezizomycetes, which includes species that typically produce large, fleshy, stipitate fruiting bodies. These fruiting bodies have rib-shaped, normally a honeycomb-shaped, bottle cap and is popular as an edible fungus. Morels commonly grow in a wide variety of habitats and are morphologically variable (Wipf et al., 1996; Dalgleish and Jacobson, 2005). Due to their unique flavor and nutritional value, these morels are used in soups and gravies, as a source of medicinal adaptogens, immunostimulants, and potential antitumor agents (Fu et al., 2013; Tietel and Masaphy, 2018). It is also widely used in traditional Chinese medicine to treat indigestion, excessive phlegm, and shortness of breath (Meng et al., 2010; Fu et al., 2013).
To date, differentially expressed genes (DEGs) between the mycelium and fruiting bodies were identified by high-throughput sequencing in several mushrooms, including Agrocybe aegerita, Auricularia polytricha, Cordyceps militaris, Ganoderma, and Lentinula edodes, and the morphological changes during the fruiting body formation were regulated by the transcription levels of development-specific genes (Song et al., 2018). Recently, the whole-genome sequence analysis of Morchella importuna and Morchella sextelata was carried out (Liu et al., 2018; Murat et al., 2018; Mei et al., 2019), but the genome and ecology of Morchella are not well understood. Besides, a series of studies indicate that the bacterial communities are associated with morel development. For example, Pseudomonas has been found to be a beneficial associate of morels previously (Pion et al., 2013); Pedobacter, Pseudomonas, Stenotrophomonas, and Flavobacterium were found to comprise the core microbiome of soil for M. sextelata ascocarps (Benucci et al., 2019). However, most bacterial community experiments were performed on the soil sample, not fruiting body samples. In our study, we conducted genome and transcriptome analysis of Morchella conica SH, as well as the endophytic bacterial communities, which will facilitate a better understanding of the molecular biology of M. conica SH.
Materials and Methods
Strains and Culture Conditions
Natural M. conica SH samples were collected from a local farm affiliated with Yunnan Academy of Agricultural Sciences, Yunnan Province, China. Fruiting bodies from different development stages (42–84 days) were randomly selected. After thorough washing with double distilled water, the fruiting bodies were immediately stored in sterilized valve bags and flash frozen in liquid nitrogen. Mycelia isolated from fruit bodies were cultivated and maintained on YPD agar (0.2% yeast extract, 0.2% peptone, 2% dextrose, and 2% agar) at 25°C in the dark.
Genome Sequencing and Assembly
High-quality genomic DNA (gDNA) was extracted from 8-day mycelia of M. conica SH using a modification of the cetyltrimethylammonium bromide (CTAB) method established by Bao et al. (2013). Optimal quality and quantity of the extracted gDNA were verified using NanoDrop (Thermo Scientific, Wilmington, DE, United States) and Qubit 2.0 fluorometer (Invitrogen, Grand Island, United States). The gDNA library for sequencing was then prepared with TruSeq DNA sample Prep Kit (Illumina, United States) according to the manufacturer’s instructions, and the library quality was performed on Agilent 2100 Bioanalyzer (Agilent, United States) and Quantiflour-ST fluorometer (Promega, United States). Whole-genome shotgun sequencing of M. conica SH was performed using an Illumina MiSeq platform (Personalbio Co., Shanghai, China). Reads with low quality (> 50% bases with Q-value ≤ 8) and Ns > 10% in read length and adapter contaminations were removed from raw data to obtain clean reads. Clean reads were then assembled using Newbler version 2.8. Gap filling was done by GapCloser v1.12 module from SOAPdenovo2 (Luo et al., 2012). The integrity of assembly was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO; Simao et al., 2015). Mitochondrial assembly was screened out from the assembly result by comparing with mitochondrial DNA of other species in Morchella using BLASTX. The Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession VOVZ00000000. The version described in this paper is VOVZ01000000.
Gene Prediction and Functional Annotation
RepeatMasker (open) v4.061 and the RepBase library2 (Bao et al., 2015) were used to mask genomic assembly scaffolds. The tRNAs were predicted by tRNAscan-SE (Lowe and Eddy, 1997) with eukaryote parameters. The rRNA sequences were aligned using BLASTN with E-value 1e–5 to identify the rRNA fragments. The snRNA genes were predicted by Erpin v5.4 software (Gautheret and Lambert, 2001) based on the Rfam database.
Genes were predicted using Augustus (Stanke et al., 2008) trained on core eukaryotic genes (CEGs) identified by Core Eukaryotic Genes Mapping Approach (CEGMA) (Parra et al., 2007). Gene model predictions were compiled using EVidenceModeler (EVM; Haas et al., 2008) including evidence from sequenced transcripts from M. conica SH assembled into contigs using PASA (Haas et al., 2011). A summary of reported Morchella genomes is listed in Table 1. Multiple sequence alignments between the M. conica SH and M. importuna CCBAS932 (JGI project ID: 1023999) genomes were performed using Nucmer (Kurtz et al., 2004) and visualized by Circos (Krzywinski et al., 2009).
Table 1. Statistics of the genome assembly of Morchella.
Annotation of predicted gene models was performed by similarity searches using BLASTp (E-value: 1e–5) based on Swiss-Prot, GenBank’s database of nonredundant proteins, and TrEMBL database (Bairoch and Apweiler, 2000). SignalP 3.03 was used to predict putative secreted proteins. KEGG Automatic Annotation Server (KAAS) analysis was used for KEGG ortholog mapping (Moriya et al., 2007), and InterProScan (Goujon et al., 2010) was used for domain annotation. Proteins were annotated using BLAST2go (Conesa et al., 2005) on the NCBI nonredundant protein database. Predicted proteases were obtained by BLASTp based on MEROPS database (Rawlings et al., 2010) and Pfam databases, respectively. To identify secondary metabolite gene clusters, the genome sequence was analyzed by online program natural product search engine (Li et al., 2009) and antiSMASH (Blin et al., 2013).
Orthology and Phylogenetic Analysis
M. conica SH and the other 15 fungi, Magnaporthe oryzae (Choi et al., 2015), Neurospora crassa FGSC 73 (Baker et al., 2015), Nectria haematococca (Coleman et al., 2009), Sclerotinia sclerotiorum (Amselem et al., 2011), Botrytis cinereal BcDW1 (Blanco-Ulate et al., 2013), Stagonospora nudorum (Hane et al., 2007), Aspergillus niger ATCC 1015 (Andersen et al., 2011), Arthrobotrys oligospora (Yang et al., 2011), Pyronema confluens (Traeger et al., 2013), Tuber melanosporum Vittad (Martin et al., 2010), Ganoderma lucidum (Chen et al., 2012), Saccharomyces cerevisiae YJM789 (Wei et al., 2007), Saccharomyces pombe (Rhind et al., 2011), M. sextelata (Mei et al., 2019), and M. importuna (Liu et al., 2018), were chosen for orthology analysis. All protein-coding genes were analyzed by OrthoMCL (Li et al., 2003) using an all-versus-all BLASTp. Genes share more than 50% sequence identity, and 50% covery was assigned into clusters. For putative orthologs or paralogs, the p-value cutoff was 1e–5. The Markov cluster (MCL) algorithm was used to assign proteins into families (Enright et al., 2002). Then, single-copy orthologs were selected to align using MUSCLE (Edgar, 2004), trimmed with Gblocks (Castresana, 2000). The maximum likelihood analysis was built by RAxML (Stamatakis, 2006) under JTT model and 100 bootstrap replicates (Burleigh et al., 2006). The divergence time between species was estimated by the Langley–Fitch method with r8s (Sanderson, 2002) calibrated based on the origin of Ascomycota at 500–600 million years ago (MYA; Taylor and Berbee, 2006).
Collinearity Analysis
All-versus-all BLASTp was performed to identify paralogous or orthologous gene pairs, and the blast result associated with M. conica SH annotation file was processed with MCScanX to identify collinear blocks (Wang et al., 2012). The results were visualized using Circos plot (Krzywinski et al., 2009).
Other Protein Families
Putative enzymes involved in carbohydrate utilization were predicted by a combination of BLASTp and HMMER (Mistry et al., 2013) searches based on carbohydrate-active enzyme database (Cantarel et al., 2009). Lipases were predicted by using BLASTp based on the Lipase Engineering Database (Fischer and Pleiss, 2003). Additionally, transporters, protein kinases, transcription factors (TFs), and lignocellulose-active proteins were classified by BLASTp (E-value cutoff 1E–10) based on Transporters Classification Database4 (Saier et al., 2006), KinBase5, Fungal Transcription Factor Database (FTFD6; Park et al., 2008), and Characterized Lignocellulose-Active Proteins of Fungal Origin Database (mycoCLAP7; Strasser et al., 2015), respectively.
Transcriptome Analysis
Gene expression profiles were evaluated in different development stages of M. conica SH. X0 is the free-living mycelium (FLM) stage, and X1, X2, X3, and X4 samples represent 1, 2, 4, and 6 weeks of fruiting body, respectively. Three biological replicates for each sample were analyzed.
Total RNA was extracted from the whole fruiting body using RNeasy Mini Kit (Qiagen, United States), according to the manufacturer’s protocols. Then, ND-2000 spectrophotometer (NanoDrop Technologies) and Bioanalyzer 2100 (Agilent Technology, United States) were used to evaluate the quantity and quality of RNA. RNA libraries were constructed and size-selected using VAHTS mRNA-seq v2 Library Prep Kit for Illumina (Vazyme, Nanjing, China) and AMPure XP Beads (Beckman, Germany), respectively. The libraries were amplified by PCR (15 cycles) and sequenced on HiSeq X platform at Shanghai Sangon Biotechnology Corporation.
Raw sequencing reads of Fastq files were filtered using seqtk8, which masked bases with Q-scores < 20. The clean data were obtained by removing reads containing adapters, ribosome RNA reads, and low-quality reads from raw data and then mapped to reference genome using Hisat2 version 2.0.4 (Kim et al., 2015). Those uniquely mapped reads were retained for read counting against the annotated genes using StringTie version 1.3.0 (Pertea et al., 2015). The expression level of the unique gene was calculated as fragments per kilobase of exon model per million fragments mapped (FPKM; Mortazavi et al., 2008) by RSEM tool (Li and Dewey, 2011). The edgeR was then used to identify DEGs (Robinson et al., 2010).
Co-expression networks were constructed using the pairwise Pearson correlations between DEGs. Louvain method for community detection was employed to part the network into modules of genes displaying similar expression profiles. To test associations between microbial taxa and modules, we used the MaAsLin2 (Morgan et al., 2012) to assess the relationship between the top 10 most abundant Morchella-associated microbes and co-expression network modules. Co-expression networks were visualized in Cytoscape (Shannon et al., 2003).
Amplicon Sample Preparation for 16S rRNA
Microbial community gDNA was extracted from different development stages of M. conica SH fruiting body using the FastDNA® Spin Kit for Soil (MP Biomedicals, GA, United States) according to manufacturer’s instructions. The DNA extract was checked on 1% agarose gel, and DNA concentration and purity were determined with NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, DE, United States). The hypervariable region V3–V4 of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) by an ABI GeneAmp® 9700 PCR thermocycler (ABI, CA, United States). The PCR amplification of 16S rRNA gene was performed as follows: initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 45 s, and single extension at 72°C for 10 min, and end at 4°C. The PCR mixtures contain 5 × TransStart FastPfu buffer 4 μl, 2.5 mM dNTPs 2 μl, forward primer (5 μM) 0.8 μl, reverse primer (5 μM) 0.8 μl, TransStart FastPfu DNA Polymerase 0.4 μl, template DNA 10 ng, and finally ddH2O up to 20 μl. PCR reactions were performed in triplicate. The PCR product was extracted from 2% agarose gel and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, United States) according to the manufacturer’s instructions and quantified using QuantusTM Fluorometer (Promega, United States).
16S rRNA Gene Sequencing and Bioinformatic Analysis
Purified amplicons were pooled in equimolar and paired-end sequenced on an Illumina MiSeq PE300 platform (Illumina, San Diego, CA, United States) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw 16S rRNA gene sequencing reads were demultiplexed, quality-filtered by fastp version 0.20.0, and merged by FLASH version 1.2.7 with the following criteria: (i) the 300-bp reads were truncated at any site receiving an average quality score of < 20 over a 50-bp sliding window, truncated reads shorter than 50 bp were discarded, and reads containing ambiguous characters were also discarded; (ii) only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence. The maximum mismatch ratio of overlap region is 0.2. Reads that could not be assembled were discarded; (iii) Samples were distinguished according to the barcode and primers, and the sequence direction was adjusted according to exact barcode matching and two nucleotides mismatch in primer matching.
Operational taxonomic units (OTUs) with 97% similarity cutoff were clustered using UPARSE version 7.1, and chimeric sequences were identified and removed. The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 based on the 16S rRNA database using a confidence threshold of 0.7.
Quantitative Real-Time PCR Analysis
Quantitative real-time PCR (qRT-PCR) was used to validate the assembled transcriptome from the RNA sequencing (RNA-seq) experiment. Total RNA was extracted from the mycelium (X0) and the whole fruiting body (X1, X2, X3, X4). RNA was isolated by using an RNA Purification Kit (TIANGen Biotech, Beijing, China). Then, total RNA (2 μg) was reverse-transcribed to cDNA by reverse transcription mix (Promega). qRT-PCR reactions were performed on an ABI 7500 system (ABI, CA, United States). The primers used in qRT-PCR are summarized in Supplementary Table 24. All qRT-PCR experiments were performed in triplicate using independent samples (Meng et al., 2019). The expression levels were determined by the 2–ΔΔCt method using CYC3 gene as internal control genes for normalization (Zhang et al., 2018).
Availability of Data
Raw transcriptomic data have been deposited in the NCBI with the SRA project ID PRJNA733830.
Results
Genome Sequencing, Assembly, and Genomic Features
The genome sequence of M. conica SH was carried out by the whole-genome shotgun sequencing strategy. A total of 7.98 Gb raw data were generated for paired-end reads. After quality filtering, 6.38-Gb clean data with 113.97 × coverage obtained (Supplementary Table 1), and 51.7-Mb genome assembly with 1.49% repeat content was ultimately used for further analysis (Supplementary Table 2). Based on K-mer statistics, the M. conica SH genome was estimated to be 55.97 Mb (Supplementary Table 3), indicating that 92.4% of the genome was captured in scaffolds. The assembly comprises 2,692 contigs (N50 = 41.7 kb) and 207 scaffolds (N50 = 0.75 Mb). According to BUSCO (Simao et al., 2015), 98.75% of the fungal single-copy orthologs were found in gene repertoires, suggesting that the assembly of M. conica SH is highly accurate. In total, 9,676 protein-coding genes were predicted, with an average transcript length of 2,336 bp and an average exon number of 4 (Table 1). The functions of these genes were predicted by different databases. Compared with other databases, large proportion genes (79.7%) were enriched in NR (Supplementary Table 4). We also found that our data are comparable to the assembly statistics in the M. importuna strains (Table 1). Furthermore, we assembled the mitochondrial genome into one scaffold with a length of 262,110 bp (Supplementary Result 1). The majority (72%) of these genes are hypothetical proteins with unknown functions. However, 590 secreted proteins (6%) were found in M. conica SH; the rest are glycoside hydrolases and proteins involved in growth and utilization of carbohydrates.
Phylogenetic Analysis of Evolutionary Relationship Between M. conica SH and Other Related Fungal Species
In order to identify the relationship between M. conica SH and other related fungal species, we compared the genome structures of the 207 scaffolds in M. conica SH assembly with that of the 504 scaffolds in the M. importuna CCBAS932 assembly (Murat et al., 2018; Figure 1A). A total of 118 matching scaffold block pairs ranging from 102 Kb to 1.4 Mb were identified between two genomes. Among those blocks, 31 contained inversions and translocations, 23 contained only inversions, 10 were completely synthetic with no inversion or translocation, and 38 contained only translocations. When comparing the predicted proteins, 7,713 of M. conica SH showed significant similarities (> 79%) to known proteins from M. snyderi DOB2414, M. importuna SCYDJ1-A1, and M. importuna CCBAS932 (Liu et al., 2018; Murat et al., 2018). Additionally, we found that most putative homologs were identified in the Pezizomycotina, and there is high sequence similarity between 4,630 predicted proteins and that from T. melanosporum (Martin et al., 2010).
Figure 1. Evolutionary relationship between Morchella conica SH and other related fungal species. (A) Genome synteny in Morchella. Circos plot showing syntenic regions in two Morchella species (M. conica SH in colored blocks, Morchella importuna in gray blocks). One syntenic block is considered only if at least two orthologous genes are present in the same order. The yellow, green, and blue links correspond to syntenic blocks composed of > 10, 5–10, and < 5 genes, respectively. (B) The phylogenetic tree of Morchella and 13 other species in the context of fungi.
To investigate the evolutionary status of M. conica SH, multigene families were predicted from proteins in M. conica SH and 14 other related representative species using the MCL algorithm (Enright et al., 2002). Finally, 672 single-copy orthologous genes were obtained and used to construct a phylogenetic tree (Supplementary Table 5). The phylogeny was inferred using RAxML version 7.1.0 with the PROTGAMMAJTT model and 100 bootstrap replicates. The divergence time of the Morchella was estimated at 25–32 MYA using r8s (Sanderson, 2002) calibrated with divergence times of 500–600 MYA for Ascomycota. M. conica and M. importuna were gathered into one cluster and separated from M. sextelata. As the sister taxon of T. melanosporum (diverged about 25–32 MYA), M. conica was significantly diverged from P. confluens (57–65 MYA), which belonged to Pezizomycetes (Figure 1B).
Functional Characterization of M. conica SH
In order to characterize the function of predicted genes, we compared the genes involved in specific pathways among M. conica SH and nine other Pezizomycetes, including five saprotrophs (M. importuna, Ascobolus immersus, Asscodesmis nigricans, Plectania confluens, and Sarcoscypha coccinea) and four ectomycorrhizal (ECM) fungi (Tuber borchii, Tuber melanosporum, Terfezia boudieri, and Choiromyces venosus), as well as 12 additional taxonomically and ecologically distinct fungi, including other ECM, saprotrophs, and pathogens (Supplementary Table 6). Genes belonging to four functional clusters were annotated.
Repertoire of Degrading Enzymes
Lignin is an aromatic macromolecule that protects the plant cellulose and hemicellulose against biodegradation. To reveal the carbohydrate-digesting capabilities of M. conica SH, we analyzed genes involved in the degradation of polysaccharides, proteins, and lipids (Figure 2 and Supplementary Tables 7–13). Overall, 232 putative proteases were identified in M. conica SH genome, the majority (171) were also found in M. importuna CCBAS932 according to the MEROPS database9 (Supplementary Table 7). Among them, 18% of proteases identified in M. conica SH contained signal peptides, the same as the pathogenic fungi (Supplementary Table 8). Additionally, P. melastoma had more proteases and secreted proteases than other species (Supplementary Tables 7, 8), while A. bisporus had the least. Compared with other members in Pezizomycetes, aspartic peptidases and threonine peptidases were significantly enriched (t-test, p < 0.005) in Morchella genomes. Notably, Pezizomycetes genomes include nearly twice the median number of lipase (Supplementary Table 9). In the M. conica SH genome, a total of 90 lipases from microsomal hydrolases family (GX-abH09) were annotated, which is significantly more (t-test, p < 0.005) than the other species. However, no difference was observed in the profile of lipases and secreted lipases among Pezizomycetes (Supplementary Tables 9, 10).
Figure 2. Double hierarchical clustering of the carbohydrate-active enzyme (CAZyme)-coding gene numbers in representative fungal genomes. A double hierarchical clustering of the number of CAZyme-coding genes for each of the fungal species was performed using the MatLab software. The Euclidian distance between gene counts was used as distance metric, and a complete linkage clustering was performed. The relative abundance of genes is represented by a color scale (on the left) from the minimum (white) to the maximum (red) number of copies per species. See Supplementary Table 6 for full names of species and lifestyles.
Ectendomycorrhiza and Ectomycorrhiza (ECM) fungal genomes contained the lowest number of carbohydrate-active enzyme (CAZyme) genes, while saprotrophic fungi contain more CAZyme genes (Supplementary Tables 11, 12). Specifically, Basidiomycota harbors more redox enzymes acting in conjunction with CAZymes (AAs) and carbohydrate esterases (CEs), while Morchella has more polysaccharide lyases (PLs) than any other fungi. In addition, saprotrophic species have nearly 50% more CAZyme genes compared with other Pezizomycetes, which possess a different lifestyle (Supplementary Figure 1A). The profiles of CAZymes among Morchella species are almost identical. CAZyme profiles revealed that M. conica SH had 177 glycoside hydrolases (GHs), 21 PLs, and 26 CEs (Supplementary Table 1A). Additionally, one cuprum peroxidase was found in the genome.
CAZyme families involved in lignocellulose degradation were analyzed in our study. Compared with ECM species, saprotroph genomes encoded a greater set of genes coding for lignocellulose oxidoreductases (Supplementary Figure 1B and Supplementary Table 13), such as cellobiose dehydrogenases (AA3), iron reductases (AA8), and lytic polysaccharide monooxygenases involving cleavage of cellulose (AA9), degradation of cellulose (GH 5_7, GH5_16, GH 6 and GH 7, GH 45 and GH 115), hemicellulose and/or pectin (GH 27, GH 28, GH 39, GH 43, GH 51, GH 53, GH 88, GH 93, PL 1, and PL 3), also glucans (GH 15 and GH 71). Among Pezizomycetes, according to BLAST analysis in lignocellulose-active protein database (see text footnote 7), 18 CAZyme families involving lignocellulose degradation were significantly enriched in Morchella (Supplementary Table 13). They are involved in the degradation of cellulose (GH 115), pectin or/and hemicellulose (GH 6, GH 26, GH 27, GH 30, GH 43, GH 79, and PL 1), and glycan (GH 12, GH 13, and GH 43).
Genes Associated With Secondary Metabolism
Plant-associated fungi produce secondary metabolites that are ecologically and nutritionally benefit to plants (Spatafora and Bushley, 2015). The core genes associated with putative secondary metabolite biosynthesis are often in gene clusters, including polyketide synthases (PKSs), non-ribosomal peptide synthetases (NRPSs), members of the 4-dimethylallyltryptophan synthetase (DMATS), and terpene synthases (TSs). Most Pezizomycetes do not contain DMATS genes compared with other types of fungi. However, P. confluens contains a large amount of secondary metabolism genes (Supplementary Table 14).
There are seven NRPS gene clusters, two PKS-like gene clusters, and one TS gene cluster in the M. conica SH genome. One of the NRPS proteins, Scaffold26.t84, has a typical extracellular siderophore domain (Supplementary Figure 2A), while another NRPS protein (Scaffold11.t41) harbors an adenylation domain (A_NRPS; cd05930), and four NRPS proteins (Scaffold7.85, Scaffold16.t118, Scaffold10.t106, and Scaffold24.t128) contain Acyl-CoA synthetase domain (CaiC; COG0318; Supplementary Figure 2B). The type I PKS gene (scaffold7.86) was adjacent to another putative NRPS gene (scaffold7.85), and the Acyl transferase domain (PksD; COG3321) was found in its product (Supplementary Figure 2C). The remaining NRPS genes have no homology to NRPS with known functions.
Transporters
In our study, a total of 269 transporter families detected in all the fungi were analyzed. As observed, 162 families were commonly presented in all species (Supplementary Table 15). Bioinformatic prediction was performed by TransportDB 2.0 (Elbourne et al., 2017), and transporter genes were clustered according to their types or functions (Supplementary Table 15). As expected, the most commonly enriched family was the major facilitator superfamily (MFS), followed by the ATP-binding cassette (ABC) and mitochondrial carrier (MC). On the other hand, BLAST analysis for the transporter classification database (TCDB; Saier et al., 2006) showed that 1,213 gene products in M. conica SH genome have homologs with > 30% identity in TCDB (Supplementary Table 16), and the most abundant is the MFS Superfamily (2.A.1). Compared with other Pezizomycetes, the membrane attack complex/perforin family (MACPF, 1.C.39) and the ion-translocating microbial rhodopsin family (MR, 3.E.1) are missing in Morchella. We also found that some transporter families are strain-specific, for example, the profiles of SH strain and CCBAS932 strain are obviously different (Supplementary Table 16). Interestingly, two transporter families, the sorting nexin27 (SNX27)-retromer assembly apparatus for recycling integral membrane proteins (SNX27-RetromerAA, 9.A.3) and the ankyrin family (8.A.28), were expanded, and the expansion is also found in C. venosus genome. However, little is known about its function in fungi. The putative ductin channel family transporter (9.A.16) was only encoded by Morchella genomes, and an adiponectin family transporter (8.A.94) was unique in M. conica SH genome. In addition, BLAST result showed that the collagen-like protein family (PFAM01391, identity > 50%) has homology in M. conica SH genome, which is related to cell adhesion, invasion, and intracellular signaling in pathogens (Humtsoe et al., 2005).
Signal Transduction
Overall, CMGC (cyclin-dependent kinase, mitogen-activated protein kinase, glycogen synthase kinase, CDC-like kinase), calcium/calmodulin-dependent kinase (CAMK), and STE (homologs of the yeast STE7, STE11, and STE20 genes) and AGC (protein kinase A, G, and C) seem to be the most common clades. Histidine kinase family was expanded in most Pezizomycetes, followed by CAMK family RAD53, which has a similar pattern in other fungi (Supplementary Table 17). Interestingly, kinase family CHK1 was the only family commonly seen in saprotrophic fungi rather than ECM fungi (t-test, p < 0.05). Among Pezizomycetes, ECM fungus T. melanosporum had the largest number of kinases (193 kinases). Several kinase families were abundant in Morchella species, such as AGC, CAMK, and STE families (t-test, p < 0.05). There are 137 protein kinases and 27 other atypical kinase genes in M. conica SH, while 153 and 25 were detected in M. importuna CCBAS932 genome, respectively. M. conica SH has 12 histidine kinases compared with 5–14 in other fungi (Supplementary Table 17).
Gene expression is controlled by the activation of different TFs, which are essential regulatory proteins involved in a variety of cell functions. Putative TFs were identified by BLASTP analysis based on the FTFD (see text footnote 6). The candidates were classified according to Park et al. (2008) (Supplementary Table 18). In general, Pezizomycetes have more TFs than other species, especially in zinc finger classes. However, there is no significant difference in the distribution of different classes within Pezizomycetes (Supplementary Figure 3). In M. conica SH, 836 TFs were predicted, including all major classes of fungal TFs, Zn2Cys6, Helix-turn-helix, AraC type, Winged helix repressor DNA-binding, and C2H2 zinc finger. Eighty-nine predicted M. conica SH proteins were similar to the characteristic TFs of T. melanosporum (BLASTp, E-value < 1e–50). Moreover, according to different cell functions, putative Morchella TFs were divided into five groups (Supplementary Table 19).
Additionally, the single mating-type locus in unifactorial Agaricomycetes was also found in the M. conica SH genome (scaffold86.t1). Individuals of heterologous fungi usually contain a MAT1 idiomorph (Turgeon and Yoder, 2000). In this study, the MAT1 gene in M. conica SH genome was homologous with MAT1-1(EAA35088) gene of N. crassa.
Comparative Transcriptome Analysis
The M. conica SH strain of the free-living mycelium (stage X0) and four different development stages of the fruiting body (X1–X4) were used for transcriptome analysis (Figure 3A and Supplementary Table 20). In total, 7,555 of 9,203 predicted genes were expressed in M. conica SH. As shown, 2,218 genes were stably expressed at each stage, 1,329 genes were specifically expressed at X0 stage, and 2,206 genes were commonly expressed in the fruiting body (Figure 3B). We then compared the genes expressed in the mycelium and fruiting body and found that 2,233 genes were expressed in both stages, while 1,329 and 2,796 genes were specifically expressed in the mycelium and fruiting body, respectively. The Gene Ontology (GO) enrichment showed that RNA polymerase II transcription activity-related genes were enriched in the mycelium-specific gene cluster, and nucleotide binding-related genes were enriched in the fruiting body-specific gene cluster (Supplementary Figure 4).
Figure 3. Expression analysis of developmental stages of Morchella conica SH. (A) Gene expression map of five developmental stages. Green and red arrows denote the numbers of significantly up- and downregulated genes [p < 0.05, fold change (FC) > 4], respectively. The development of X0 (mycelium), X1–X2 young fruiting body (YFB), and X3–X4 fruiting body (FB) is shown. Each sample was sequenced in three biological replicates. (B) Venn plot shows the number of detected genes overlapping in the compared developmental stages. (C) Heatmap shows the differentially expressed genes (DEGs) during the five development stages of M. conica SH; group 1 and group 2 represent distinct expression patterns. (D). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of genes in group 1 and group 2.
The DEGs (adjusted p < 0.05, fold change ≥ 2) between different development stages were identified. Comparing between X0–X1, X1–X2, X2–X3, and X3–X4 growth stages, 1,125, 396, 551, and 650 protein-coding genes were significantly up- or down-regulated (p < 0.05, fold change ≥ 2), respectively (Figure 3 and Supplementary Table 20). To study the function of the DEGs in the five development stages, GO terms and KEGG Orthology (KO) terms were used to classify genes to functional categories and pathways. DEGs were grouped into 16 clusters by parasitic dynamics using the K-Means clustering algorithm (Figure 4). Genes with a high expression level in the mycelium were enriched in clusters 4, 5, and 7, and GO analysis displays that they were involved in cellular component and transport, while genes with high expression levels in the fruiting body were enriched in clusters 2, 3, 10, and 13 and related to phospholipid metabolic process and kinase activity (Supplementary Table 21). The expression pattern was illustrated by Heatmap. Two groups of genes showed distinct expression patterns (Figure 3D). Genes belonging to group 1 were significantly upregulated in the fruiting body (X1–X4), and KEGG analysis suggested that they were mainly involved in other glycan degradation and ABC transporters, while genes belonging to group 2 were dramatically downregulated in the fruiting body and were mainly enriched in taurine and hypotaurine metabolism, tyrosine metabolism, and styrene degradation.
Figure 4. Dynamic progression of Morchella conica SH transcriptome in different development stages. K-means clustering showing the expression profile of the M. conica SH transcriptome. Sixteen clusters were identified and presented along the five stages (X0, X1, X2, X3, and X4) from 9,203 differentially expressed genes; genes of each cluster were listed in Supplementary Table 21.
In order to identify core genes involved in biological processes during the development of the fruiting body, we compared transcriptome profiles of undifferentiated mycelium grown on agar medium with the fruiting body at different stages. Of the 7,555 transcripts detected in the fruiting body, 226 (3%) were significantly upregulated (p < 0.05) in all the stages of the fruit body. Among the 50 most highly induced genes in the fruiting body, transporters, fungal transcriptional factors, and kinases were among the most highly induced genes (Supplementary Table 22). Several transcripts coding for proteases (M24A, T01A, I87, M18, and T01X families) were upregulated transcripts during X0–X1, and the most upregulated protease was one of T01A family (Protein ID Scaffold23.t27; Supplementary Table 23), showing a 300-fold induction in the fruiting body. Ten proteases (C26, S08A, M24X, T03, A01A, C12, M28X, and M41 families) were downregulated in the fruiting body, including the gene encoding deubiquitinating enzyme Uch2 (scaffold2.t334) from C12 family. Among the lipases, an Alpha esterase gene (scaffold15.t73) from abH01 family was significantly upregulated in the fruiting body, while scaffold47.t53 and scaffold39.t45 showed a 9,000-fold and 100-fold induction in the mycelium.
Validation of Gene Expression by qRT-PCR
To validate the RNA-Seq results, six genes were randomly selected from the significant gene expression patterns of group 1 and group 2 and analyzed by qRT-PCR. The CYC3 gene was used as an internal reference gene for qRT-PCR analyses, which has been evaluated as the most stable gene in different development stages (Zhang et al., 2018). For each gene, the expression count values of transcriptome data exhibited a similar expression profile compared with the qPCR results, and the RNA-seq and qRT-PCR results revealed a strong correlation with each other (Pearson correlation, r = 0.8715, p = 3.5421 E- 15), suggesting reliable expression results generated via RNA-seq (Supplementary Figure 6).
Dynamic Shift of Microbial Communities in the Development of the Fruiting Body
Microbial communities can trigger the formation of primordia and the development of the fruiting body. To understand the microbial communities in M. conica SH, we examined the microbial communities in four development stages of the fruiting body. By detecting bacterial 16S rRNA, 515 OTUs were identified in at least one stage (Supplementary Table 25). When comparing different stages, we found that no distinct changes on the Chao index of OTU levels during X1–X3, but it was significantly increased at X4 (Figure 5A). Venn plot showed that a large number of bacterial communities (106) can be detected in the four stages, which might be the essential bacteria for the development of the fruiting body. And 18, 26, 18, and 152 specific communities were detected at the four stages, suggesting more communities might be involved in the degradation of useless substrate at the mature stage (Figure 5B). To further identify the dominant bacterial communities at each stage, we calculated the percentage of community abundance on the genus level. As observed, high abundances of Flavobacterium, Pseudomonas, and Pedobacter were found in the four stages, while Massilia and Herbaspirillum were highly abundant at the first three stages or mature stages, respectively (Figure 5C). As previously reported, Flavobacterium associated with promoting the development of its host (Chaudhary et al., 2019; Kim and Yu, 2020). Massilia was involved in the synthesis of multiple secondary metabolites and enzymes, phosphorus solubilization, degradation of phenanthrene, and resistance to heavy metals (Feng et al., 2016; Zheng et al., 2017). Herbaspirillum could fix nitrogen in rice and sugarcane (Hoseinzade et al., 2016). Our results suggested that these bacterial communities might play an important role in the development and maturity of M. conica SH.
Figure 5. Different microbial communities in the fruiting body development stages. (A) The Chao index of operational taxonomic unit (OTU) levels. (B) Venn plot shows the number of detected OTUs of the fruiting body. (C) The endophytic microbial community analysis with development stages of the fruiting body.
In order to check the possible regulatory roles of bacterial communities in gene expression, we investigated the GO functions of bacteria-regulated DEGs. Three main gene clusters highly correlated with the abundance of bacterial communities were found in our analysis. GO enrichment displayed genes associated with Arsenicitalea and Duganella involved in fungal-type cell wall, cell wall, external encapsulating structure, cell surface, and hyphal cell wall; genes associated with Herbaspirillum, Pedobacter, Devosia, and Stenotrophomonas involved in hyphal cell wall, proteasome complex, ubiquitin ligase complex, nuclear ubiquitin ligase complex, and TF TFIIE complex; and genes associated with Chryseobacterium involved in proteasome regulatory particle, proteasome accessory complex, and integral and intrinsic component of plasma membrane (Figure 6). Taken together, our results suggested that endophytic bacterial communities might promote the development of the fruiting body through regulating specific gene expression in M. conica SH.
Figure 6. Gene Ontology (GO) functions of bacteria-regulated differentially expressed genes (DEGs). Using fruiting body-associated DEGs (gray dots), we created a co-expression network, which was partitioned into modules by network community detection (colored lines enclosed). Top five enriched GO terms were indicated for each network module. Hypergeometric tests revealed the bacteria-regulated genes (colored dots) that mapped the network.
Discussion
M. conica was not only consumed as tasteful food, health nutritional supplement for its high gastronomic quality, fatigue resistance, and gastroprotective effects but also used in traditional Chinese medicine to potentially treat indigestion, excessive phlegm, shortness of breath, and antitumor and immunomodulatory activities (Tietel and Masaphy, 2018). However, the understanding of its basic biology and genome information is very limited. In this study, we reported the complete genome blueprint of M. conica SH and constructed a transcriptome and endophytic bacterial community analysis of M. conica SH.
This study is the first genome sequencing analysis of M. conica SH, whose general assembly features are similar to those of other reported Morchella species (Table 1). Comparison between M. conica SH and M. importuna CCBAS932 in genome structures showed inversions and translocations occurring in several homologous scaffolds, suggesting that evolutionary genome rearrangement accrued between two species. The ecology of Morchella spp. is not well understood; some species showed a symbiotic relationship with plant root, while some were saprotrophs (Murat et al., 2018). We compared the genome content of M. conica SH with those of M. importuna CCBAS932 and 15 other Pezizomycetes. Our research revealed a similar distribution of protein families between M. conica SH and other saprotrophic fungi, suggesting that M. conica has a saprotrophic lifestyle. Many saprotrophic fungi resemble pathogens in their repertoire of degrading enzymes; the secreted proteases in pathogenic fungi are mainly involved in the infection process (Lucking et al., 2009). However, in M. conica SH and other saprotrophs, they are more likely to facilitate adaption on complex cultivated environments (Heneghan et al., 2009; Erjavec et al., 2012). The most abundant protease family in M. conica SH is serine protease, which belongs to fibrinolytic enzymes in other edible fungi. They have great potential to be developed as functional food additives for the prevention and treatment of thrombotic diseases (Liu et al., 2017). Thirteen GH families were found in M. conica SH, including the GH61 family. It was reported that its cellulolytic activity was weak, and it was also detected in T. melanosporum Vittad (Karkehabadi et al., 2008; Martin et al., 2010). Ligninolytic fungi degrade lignin using combinations of multiple isoenzymes of three heme peroxidases: lignin peroxidases, manganese peroxidases, and hybrid enzymes (Zhao et al., 2013). The genome of M. conica SH contains a putative manganese peroxidase and lignin peroxidase but lacks genes encoding hybrid enzymes. Many morels (Morchella) were strongly biased toward burned tree trunks after wildfire (Greene et al., 2010); damaged trees may require less enzymes to degrade M. conica SH.
The transporter system plays a critical role in fundamental cellular processes, which are distributed through membranes to deliver essential nutrients, excrete waste products, and help cell signal transduction under environmental changes. MFS and ABC are the most common types of transporters in all sequenced fungal genomes and are usually involved in the transport of different substrates (Morschhauser, 2010). ABCs are important in defending pathogens from toxic compounds produced by the host. The profile of secondary transporter genes can reflect the fungal lifestyle. For example, Pezizomycete genomes lack several amino acid and ion transporter genes, but they encode other families to compensate (Supplementary Table 16). Potassium ions can be transported via KUP (K+ uptake) and Trk (K+ transporter) families (Schlosser et al., 1995). In ECMs, many species lack KUP family proteins, but they could use Trk family K+:H+ symporter. In contrast, saprotrophic fungi are enriched in genes encoding potassium ion transporters. Furthermore, MACPFs are restricted to certain species of Pezizomycotina and Basidiomycota and have multiple biological roles in cell adhesion and signaling (Humtsoe et al., 2005).
Based on our genome sequencing and gene annotation, we performed transcriptome analysis of M. conica SH in different development stages and mainly have three new discoveries: 1) RNA polymerase II TF activity was enriched in the mycelium specifically enriched gene cluster. As known, the initiation of pol II transcription is of particular interest because of the regulation of the process. Transcriptional activators and repressors exert their effects at this early stage in gene expression to influence cellular physiology and development (Liu et al., 2013). The enrichment of transcription system in the mycelium indicated that most genes might be active in the first stages and repressed in late stages rather than active in late stages. 2) The transition from vegetative mycelium to mature fruit body has undergone significant morphological changes. Genes encoding hydrophobins, CAZymes, and transcriptional factors were proven to be essential in this biological progress. In Morchella, genes involved in the vegetative growth of M. importuna have been investigated previously and DEGs have been found in primary metabolism processes, such as carbohydrates, fatty acids, proteins, and energy metabolism. Major developmental switch started in an earlier stage of fruit body growth, and less differences were found in later stages (Liu et al., 2019). KEGG analysis of genes significantly upregulated in the fruiting body suggested that they were mainly involved in other glycan degradation and ABC transporters. MFS and ABC transporters are usually involved in the transport of different substrates and important in defending pathogen (Morschhauser, 2010). The upregulation of these genes suggested that increasing disease resistance plays an important role in the development. 3) We also checked the genes annotated in our study and found that the downregulated protease gene was deubiquitinating enzyme Uch2 (scaffold2.t334) from C12 family, indicating an 850-fold reduction in fruit body. Many deubiquitinating enzymes are associated with the 26S proteasome contributing to the regulation of the activity of particle, and they are critical in the response to stresses such as nutrient starvation and heat shock (Komander et al., 2009). Attributed to the low-nutrition environment during the growth of fruit body, the deubiquitinating enzyme may regulate stress responses and the corresponding signal transductions. Among the lipases, Alpha esterase gene (scaffold15.t73) from abH01 family was significantly upregulated in fruit body, and the function of those epoxide hydrolases in fungi was unknown. They may play an important role in lipid metabolism during mycelium growth and providing M. conica SH the capability of degrading its host cell walls during colonization. On the other hand, kinase Akt from AGC family was enriched in Morchella species. Specifically, family plays a critical role in cell growth and survival (Kane and Weiss, 2003), and fungi employ PKA pathway in a variety of processes including control of differentiation and sexual development (Xu and Hamer, 1996). Furthermore, when comparing the TFs between T. melanosporum and M. conica SH, only three homologous TFs (TmelHMS1, TmelMBF1, and TmelHAA1) in M. conica SH showed significant upregulation in the fruiting body (fold change > 200). As known, TmelHMS1 functions as a regulator of cell morphology and TmelHaa1 as a weak acid stress response regulator. This is in contrast to what is observed in T. melanosporum: the TmeHMS1 was upregulated in ectomycorrhizas instead of fruit body (Montanini et al., 2011), indicating this homolog may be a critical and common regulator in different morphologies. Furthermore, the formation of the fruiting body is a highly complex developmental process, and the mating-type loci are the master regulators in fungi (Montanini et al., 2011). Homologs (Scaffold86.t1) of a fruit body regulator protein MAT-1-2-1 (Martin et al., 2010) are present in the M. conica SH genome; however, the expression levels of these genes are low in both the mycelium and fruit body.
Many studies have noted beneficial interactions between bacteria and other mushrooms; endophytes exist at a certain or the entire stage of the host and can facilitate the growth and enhance the stress resistance of the host. Previous research mainly focused on the bacterial communities in the soil; in this study, we firstly detect the endophytic bacterial communities in different development stages of the fruiting body in M. conica SH. Our results showed that Flavobacterium, Pseudomonas, and Pedobacter were dominant in the whole development stages of the fruiting bodies. The high abundance of these endophytic bacterial communities raises questions concerning their roles in the development of M. conica SH. It is reported that Pedobacter, Pseudomonas, Stenotrophomonas, and Flavobacterium were dominant in the microbiome of M. sextelata cultivated soil (Benucci et al., 2019), and Pseudomonas, Flavobacterium, Janthinobacterium, and Polaromonas were also detected in the fruiting bodies of Pezizales truffle species (Benucci and Bonito, 2016; Splivallo et al., 2019), including Kalapuya brunnea, which belongs to the Morchellaceae family (Trappe et al., 2010). The role of Flavobacterium, Pseudomonas, and Pedobacter in the development of indicated species indicated that they were conserved in the Morchella. Additionally, Massilia and Herbaspirillum showed distinct patterns in the first three stages and mature stages. As reported, Massilia was involved in the synthesis of multiple secondary metabolites and enzymes (Feng et al., 2016; Zheng et al., 2017), while Herbaspirillum was involved in nitrogen fixation (Hoseinzade et al., 2016). The obvious differences in abundance and functions of these two endophytic bacterial communities provide a direction for our future research. We will investigate the amount of different endophytic bacterial communities in the early and late stages and identification of the function of endophytic bacteria. Then, by cultivating endophytic bacteria with specific functions such as promoting growth, resisting pathogens, or producing active substances (Yuan et al., 2018; Ullah et al., 2019) by putting these bacteria back (injection or co-cultivation) into the host, they could complete their special function (Afzal et al., 2019; Eke et al., 2019).
Endophytic bacteria play many important roles during the whole life of the host, with biological effects on plant growth and development and disease resistance (Ryan et al., 2008). Endophytic bacteria promote plant growth and development by benefiting plants in getting nutrients, and the other is to modulate the production of growth and development-related hormones (Ma et al., 2016). Endophytic bacteria can protect plants from the attack of pathogens by inhibiting the growth of plant pathogens, such as the production of antibiotics and lysozymes (Miliute et al., 2015), causing the pathogens to lack nutrients and initiating plant defense mechanisms (Afzal et al., 2019).
Finally, in order to further verify the role of endophytic bacterial communities in the growth and development of M. conica SH, we analyzed the relationship between the abundance of endophytic bacterial communities and gene expression at different development stages of the fruiting body and found that several endophytic bacterial communities, including Herbaspirillum, Pedobacter, Arsenicitalea, and Duganella, showed high correlation with modules enriched for “cell wall” and “proteasome regulatory subunit”; cell wall synthesis and remodeling are a common function during the fruiting body development (Krizsan et al., 2019). Proteasome is regarded as an important regulatory mechanism in cell cycle and growth and is highly expressed in several fruiting body-formatted fungi (Yamada et al., 2006; Wang et al., 2013). Taken together, our findings suggested that bacterial communities might promote the development of the fruiting body through regulating specific gene expressions in pathways relating to cell wall synthesis and proteasome of M. conica SH. Further studies are warranted for unraveling the potential role of endophytic bacteria in the fruiting body development.
Conclusion
In this study, we first reported the genome and transcriptomes of the medicinal edible fungi M. conica SH and provide essential tools to unravel the mechanisms of substrate degradation and fruit body formation. Transcriptome analysis reveals significant differences between the mycelia and fruiting body. ABC transporters and metabolism-related genes are responsible for the mycelia and fruiting body, respectively. Besides, we first detected the endophytic bacterial communities in the fruiting body. Dynamic shifts of bacterial communities during different development stages of the fruiting body were detected, and co-expression analysis suggested that bacterial communities might play an important role in regulating specific gene expressions. In summary, our study provided a direction for artificial cultivation of M. conica SH on gene expression and bacterial communities.
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 in the article/Supplementary Material.
Author Contributions
BL, GW, YS, LZ, and XiW analyzed the genomic data and were responsible for the experiment. PL, JW, HL, and AP analyzed the metagenomic information. YZ, XL, and WC provided the samples. XL and WC validated of gene expression by qRT-PCR. BL, GW, YS, and QM wrote the manuscript. LS and YH edited the manuscript. YY, CZ, and XuW analyzed the mitochondrial genome information of M. conica SH and added orthology analysis of M. importuna and M. sextelata. CZ and XuW analyzed the genome information and edited the revised manuscript. XT conceived the experimental design and edited the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the Agriculture Research System of Shanghai, China (Grant No. 201710), Shanghai Science and Technology Commission The Belt and Road Project (20310750500), Key Projects of Shanghai Science and Technology Commission (19DZ1203802), Shanghai Agricultural Science and Technology Innovation Action Plan (21N31900800), Shanghai Academy of Agricultural Sciences Youth Talent Program (ZP21211), 2021 SAAS Project on Agricultural Science and Technology Innovation Supporting Area [SAAS Application Basic Study 2021 (09)], 2021 SAAS Project on Agricultural Science and Technology Innovation Supporting Area [SAAS Platform 2021 (08)], and Project of Renovation Capacity Building for the Young Sci-Tech Talents Sponsored by Xinjiang Academy of Agricultural Sciences (xjnkq-2019016).
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.
Acknowledgments
All authors wish to thank Z. P. Zhou, J. L. Chen, and M. J. Chen for their assistance with processing and uploading data.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.682356/full#supplementary-material
Footnotes
- ^ http://www.repeatmasker.org/
- ^ https://www.girinst.org/
- ^ http://www.cbs.dtu.dk/services/SignalP/
- ^ http://www.tcdb.org
- ^ http://kinase.com/
- ^ http://ftfd.snu.ac.kr
- ^ https://mycoclap.fungalgenomics.ca
- ^ http://github.com/lh3/seqtk
- ^ http://merops.sanger.ac.uk
References
Afzal, I., Shinwari, Z. K., Sikandar, S., and Shahzad, S. (2019). Plant beneficial endophytic bacteria: mechanisms, diversity, host range and genetic determinants. Microbiol. Res. 221, 36–49. doi: 10.1016/j.micres.2019.02.001
Amselem, J., Cuomo, C. A., Van Kan, J. A., Viaud, M., Benito, E. P., Couloux, A., et al. (2011). Genomic analysis of the necrotrophic fungal pathogens Sclerotinia sclerotiorum and Botrytis cinerea. PLoS Genet. 7:e1002230. doi: 10.1371/journal.pgen.1002230
Andersen, M. R., Salazar, M. P., Schaap, P. J., Van De Vondervoort, P. J., Culley, D., Thykaer, J., et al. (2011). Comparative genomics of citric-acid-producing Aspergillus niger ATCC 1015 versus enzyme-producing CBS 513.88. Genome Res. 21, 885–897. doi: 10.1101/gr.112169.110
Bairoch, A., and Apweiler, R. (2000). The SWISS-PROT protein sequence database and its supplement TrEMBL in 2000. Nucleic Acids Res. 28, 45–48. doi: 10.1093/nar/28.1.45
Baker, S. E., Schackwitz, W., Lipzen, A., Martin, J., Haridas, S., Labutti, K., et al. (2015). Draft genome sequence of Neurospora crassa Strain FGSC 73. Genome Announc. 3:e00074-15. doi: 10.1128/genomeA.00074-15
Bao, D., Gong, M., Zheng, H., Chen, M., Zhang, L., Wang, H., et al. (2013). Sequencing and comparative analysis of the straw mushroom (Volvariella volvacea) genome. PLoS One 8:e58294. doi: 10.1371/journal.pone.0058294
Bao, W., Kojima, K. K., and Kohany, O. (2015). Repbase Update, a database of repetitive elements in eukaryotic genomes. Mob. DNA 6:11. doi: 10.1186/s13100-015-0041-9
Benucci, G. M. N., and Bonito, G. M. (2016). The Truffle Microbiome: Species and Geography Effects on Bacteria Associated with Fruiting Bodies of Hypogeous Pezizales. Microb Ecol. 72, 4–8. doi: 10.1007/s00248-016-0755-3
Benucci, G. M. N., Longley, R., Zhang, P., Zhao, Q., Bonito, G., and Yu, F. Q. (2019). Microbial communities associated with the black morel Morchella sextelata cultivated in greenhouses. PeerJ 7:19. doi: 10.7717/peerj.7744
Blanco-Ulate, B., Allen, G., Powell, A. L., and Cantu, D. (2013). Draft genome sequence of Botrytis cinerea BcDW1, inoculum for noble rot of grape berries. Genome Announc. 1:e00252-13. doi: 10.1128/genomeA.00252-13
Blin, K., Medema, M. H., Kazempour, D., Fischbach, M. A., Breitling, R., Takano, E., et al. (2013). antiSMASH 2.0–a versatile platform for genome mining of secondary metabolite producers. Nucleic Acids Res. 41, W204–W212. doi: 10.1093/nar/gkt449
Burleigh, J. G., Driskell, A. C., and Sanderson, M. J. (2006). Supertree bootstrapping methods for assessing phylogenetic variation among genes in genome-scale data sets. Syst. Biol. 55, 426–440. doi: 10.1080/10635150500541722
Cantarel, B. L., Coutinho, P. M., Rancurel, C., Bernard, T., Lombard, V., and Henrissat, B. (2009). The Carbohydrate-Active EnZymes database (CAZy): an expert resource for Glycogenomics. Nucleic Acids Res. 37, D233–D238. doi: 10.1093/nar/gkn663
Castresana, J. (2000). Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552. doi: 10.1093/oxfordjournals.molbev.a026334
Chaudhary, D. K., Kim, D.-U., Kim, D., and Kim, J. (2019). Flavobacterium petrolei sp. nov., a novel psychrophilic, diesel-degrading bacterium isolated from oil-contaminated Arctic soil. Sci. Rep. 9:4134. doi: 10.1038/s41598-019-40667-7
Chen, S., Xu, J., Liu, C., Zhu, Y., Nelson, D. R., Zhou, S., et al. (2012). Genome sequence of the model medicinal mushroom Ganoderma lucidum. Nat. Commun. 3:913. doi: 10.1038/ncomms1923
Choi, J., Chung, H., Lee, G. W., Koh, S. K., Chae, S. K., and Lee, Y. H. (2015). Genome-Wide analysis of hypoxia-responsive genes in the rice blast fungus, Magnaporthe oryzae. PLoS One 10:e0134939. doi: 10.1371/journal.pone.0134939
Coleman, J. J., Rounsley, S. D., Rodriguez-Carres, M., Kuo, A., Wasmann, C. C., Grimwood, J., et al. (2009). The genome of Nectria haematococca: contribution of supernumerary chromosomes to gene expansion. PLoS Genet. 5:e1000618. doi: 10.1371/journal.pgen.1000618
Conesa, A., Gotz, S., Garcia-Gomez, J. M., Terol, J., Talon, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Dalgleish, H. J., and Jacobson, K. M. (2005). A first assessment of genetic variation among Morchella esculenta (morel) populations. J. Hered. 96, 396–403. doi: 10.1093/jhered/esi045
Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797. doi: 10.1093/nar/gkh340
Eke, P., Kumar, A., Sahu, K. P., Wakam, L. N., Sheoran, N., Ashajyothi, M., et al. (2019). Endophytic bacteria of desert cactus (Euphorbia trigonas Mill) confer drought tolerance and induce growth promotion in tomato (Solanum lycopersicum L.). Microbiol. Res. 228:126302. doi: 10.1016/j.micres.2019.126302
Elbourne, L. D., Tetu, S. G., Hassan, K. A., and Paulsen, I. T. (2017). TransportDB 2.0: a database for exploring membrane transporters in sequenced genomes from all domains of life. Nucleic Acids Res. 45, D320–D324. doi: 10.1093/nar/gkw1068
Enright, A. J., Van Dongen, S., and Ouzounis, C. A. (2002). An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res. 30, 1575–1584. doi: 10.1093/nar/30.7.1575
Erjavec, J., Kos, J., Ravnikar, M., Dreo, T., and Sabotic, J. (2012). Proteins of higher fungi–from forest to application. Trends Biotechnol. 30, 259–273. doi: 10.1016/j.tibtech.2012.01.004
Feng, G. D., Yang, S. Z., Li, H. P., and Zhu, H. H. (2016). Massilia putida sp nov., a dimethyl disulfide-producing bacterium isolated from wolfram mine tailing. Int. J. Syst. Evol. Microbiol. 66, 50–55. doi: 10.1099/ijsem.0.000670
Fischer, M., and Pleiss, J. (2003). The lipase engineering database: a navigation and analysis tool for protein families. Nucleic Acids Res. 31, 319–321. doi: 10.1093/nar/gkg015
Fu, L., Wang, Y., Wang, J., Yang, Y., and Hao, L. (2013). Evaluation of the antioxidant activity of extracellular polysaccharides from Morchella esculenta. Food Funct. 4, 871–879. doi: 10.1039/c3fo60033e
Gautheret, D., and Lambert, A. (2001). Direct RNA motif definition and identification from multiple sequence alignments using secondary structure profiles. J. Mol. Biol. 313, 1003–1011. doi: 10.1006/jmbi.2001.5102
Goujon, M., Mcwilliam, H., Li, W., Valentin, F., Squizzato, S., Paern, J., et al. (2010). A new bioinformatics analysis tools framework at EMBL-EBI. Nucleic Acids Res. 38, W695–W699. doi: 10.1093/nar/gkq313
Greene, D. F., Hesketh, M., and Pounden, E. (2010). Emergence of morel (Morchella) and pixie cup (Geopyxis carbonaria) ascocarps in response to the intensity of forest floor combustion during a wildfire. Mycologia 102, 766–773. doi: 10.3852/08-096
Haas, B. J., Salzberg, S. L., Zhu, W., Pertea, M., Allen, J. E., Orvis, J., et al. (2008). Automated eukaryotic gene structure annotation using EVidenceModeler and the Program to Assemble Spliced Alignments. Genome Biol. 9:R7. doi: 10.1186/gb-2008-9-1-r7
Haas, B. J., Zeng, Q., Pearson, M. D., Cuomo, C. A., and Wortman, J. R. (2011). Approaches to fungal genome annotation. Mycology 2, 118–141. doi: 10.1080/21501203.2011.606851
Hane, J. K., Lowe, R. G., Solomon, P. S., Tan, K. C., Schoch, C. L., Spatafora, J. W., et al. (2007). Dothideomycete plant interactions illuminated by genome sequencing and EST analysis of the wheat pathogen Stagonospora nodorum. Plant Cell 19, 3347–3368. doi: 10.1105/tpc.107.052829
Heneghan, M. N., Porta, C., Zhang, C., Burton, K. S., Challen, M. P., Bailey, A. M., et al. (2009). Characterization of serine proteinase expression in Agaricus bisporus and Coprinopsis cinerea by using green fluorescent protein and the A. bisporus SPR1 promoter. Appl. Environ. Microbiol. 75, 792–801. doi: 10.1128/AEM.01897-08
Hoseinzade, H., Ardakani, M. R., Shahdi, A., Rahmani, H. A., Noormohammadi, G., and Miransari, M. (2016). Rice (Oryza sativa L.) nutrient management using mycorrhizal fungi and endophytic Herbaspirillum seropedicae. J. Integr. Agric. 15, 1385–1394. doi: 10.1016/s2095-3119(15)61241-2
Humtsoe, J. O., Kim, J. K., Xu, Y., Keene, D. R., Hook, M., Lukomski, S., et al. (2005). A streptococcal collagen-like protein interacts with the alpha2beta1 integrin and induces intracellular signaling. J. Biol. Chem. 280, 13848–13857. doi: 10.1074/jbc.M410605200
Kane, L. P., and Weiss, A. (2003). The PI-3 kinase/Akt pathway and T cell activation: pleiotropic pathways downstream of PIP3. Immunol. Rev. 192, 7–20. doi: 10.1034/j.1600-065x.2003.00008.x
Karkehabadi, S., Hansson, H., Kim, S., Piens, K., Mitchinson, C., and Sandgren, M. (2008). The first structure of a glycoside hydrolase family 61 member, Cel61B from Hypocrea jecorina, at 1.6 A resolution. J. Mol. Biol. 383, 144–154. doi: 10.1016/j.jmb.2008.08.016
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317
Kim, H., and Yu, S. M. (2020). Flavobacterium nackdongense sp. nov., a cellulose-degrading bacterium isolated from sediment. Arch. Microbiol. 202, 591–595. doi: 10.1007/s00203-019-01770-5
Komander, D., Clague, M. J., and Urbe, S. (2009). Breaking the chains: structure and function of the deubiquitinases. Nat. Rev. Mol. Cell Biol. 10, 550–563. doi: 10.1038/nrm2731
Krizsan, K., Almasi, E., Merenyi, Z., Sahu, N., Viragh, M., Koszo, T., et al. (2019). Transcriptomic atlas of mushroom development reveals conserved genes behind complex multicellularity in fungi. Proc. Natl. Acad. Sci. U.S.A. 116, 7409–7418. doi: 10.1073/pnas.1817822116
Krzywinski, M., Schein, J., Birol, I., Connors, J., Gascoyne, R., Horsman, D., et al. (2009). Circos: an information aesthetic for comparative genomics. Genome Res. 19, 1639–1645. doi: 10.1101/gr.092759.109
Kurtz, S., Phillippy, A., Delcher, A. L., Smoot, M., Shumway, M., Antonescu, C., et al. (2004). Versatile and open software for comparing large genomes. Genome Biol. 5:R12. doi: 10.1186/gb-2004-5-2-r12
Li, B., and Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323. doi: 10.1186/1471-2105-12-323
Li, L., Stoeckert, C. J. Jr., and Roos, D. S. (2003). OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 13, 2178–2189. doi: 10.1101/gr.1224503
Li, M. H., Ung, P. M., Zajkowski, J., Garneau-Tsodikova, S., and Sherman, D. H. (2009). Automated genome mining for natural products. BMC Bioinformatics 10:185. doi: 10.1186/1471-2105-10-185
Liu, W., Cai, Y., He, P., Chen, L., and Bian, Y. (2019). Comparative transcriptomics reveals potential genes involved in the vegetative growth of Morchella importuna. 3 Biotech 9:81. doi: 10.1007/s13205-019-1614-y
Liu, W., Chen, L., Cai, Y., Zhang, Q., and Bian, Y. (2018). Opposite polarity monospore genome de novo sequencing and comparative analysis reveal the possible heterothallic life cycle of Morchella importuna. Int. J. Mol. Sci. 19:2525. doi: 10.3390/ijms19092525
Liu, X., Bushnell, D. A., and Kornberg, R. D. (2013). RNA polymerase II transcription: structure and mechanism. Biochim. Biophys. Acta Gene Regul. Mech. 1829, 2–8. doi: 10.1016/j.bbagrm.2012.09.003
Liu, X., Kopparapu, N.-K., Li, Y., Deng, Y., and Zheng, X. (2017). Biochemical characterization of a novel fibrinolytic enzyme from Cordyceps militaris. Int. J. Biol. Macromol. 94, 793–801. doi: 10.1016/j.ijbiomac.2016.09.048
Lo, Y. C., Lin, S. Y., Ulziijargal, E., Chen, S. Y., Chien, R. C., Tzou, Y. J., et al. (2012). Comparative study of contents of several bioactive components in fruiting bodies and mycelia of culinary-medicinal mushrooms. Int. J. Med. Mushrooms 14, 357–363. doi: 10.1615/intjmedmushr.v14.i4.30
Lowe, T. M., and Eddy, S. R. (1997). tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 25, 955–964. doi: 10.1093/nar/25.5.955
Lucking, R., Huhndorf, S., Pfister, D. H., Plata, E. R., and Lumbsch, H. T. (2009). Fungi evolved right on track. Mycologia 101, 810–822. doi: 10.3852/09-016
Luo, R., Liu, B., Xie, Y., Li, Z., Huang, W., Yuan, J., et al. (2012). SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience 1:18. doi: 10.1186/2047-217X-1-18
Ma, Y., Rajkumar, M., Zhang, C., and Freitas, H. (2016). Beneficial role of bacterial endophytes in heavy metal phytoremediation. J. Environ. Manage. 174, 14–25. doi: 10.1016/j.jenvman.2016.02.047
Martin, F., Kohler, A., Murat, C., Balestrini, R., Coutinho, P. M., Jaillon, O., et al. (2010). Perigord black truffle genome uncovers evolutionary origins and mechanisms of symbiosis. Nature 464, 1033–1038. doi: 10.1038/nature08867
Mei, H., Qingshan, W., and Baiyintala, and Wuhanqimuge. (2019). The whole-genome sequence analysis of Morchella sextelata. Sci. Rep. 9:15376. doi: 10.1038/s41598-019-51831-4
Meng, F., Zhou, B., Lin, R., Jia, L., Liu, X., Deng, P., et al. (2010). Extraction optimization and in vivo antioxidant activities of exopolysaccharide by Morchella esculenta SO-01. Bioresour. Technol. 101, 4564–4569. doi: 10.1016/j.biortech.2010.01.113
Meng, J., Chen, X., and Zhang, C. (2019). Transcriptome-based identification and characterization of genes responding to imidacloprid in Myzus persicae. Sci. Rep. 9:13285. doi: 10.1038/s41598-019-49922-3
Miliute, I., Buzaite, O., Baniulis, D., and Stanys, V. (2015). Bacterial endophytes in agricultural crops and their role in stress tolerance: a review. Þemdirbystë (Agriculture) 102, 465–478. doi: 10.13080/z-a.2015.102.060
Mistry, J., Finn, R. D., Eddy, S. R., Bateman, A., and Punta, M. (2013). Challenges in homology search: HMMER3 and convergent evolution of coiled-coil regions. Nucleic Acids Res. 41:e121. doi: 10.1093/nar/gkt263
Montanini, B., Levati, E., Bolchi, A., Kohler, A., Morin, E., Tisserant, E., et al. (2011). Genome-wide search and functional identification of transcription factors in the mycorrhizal fungus Tuber melanosporum. New Phytol. 189, 736–750. doi: 10.1111/j.1469-8137.2010.03525.x
Morgan, X. C., Tickle, T. L., Sokol, H., Gevers, D., Devaney, K. L., Ward, D. V., et al. (2012). Dysfunction of the intestinal microbiome in inflammatory bowel disease and treatment. Genome Biol. 13:R79. doi: 10.1186/gb-2012-13-9-r79
Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A. C., and Kanehisa, M. (2007). KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 35, W182–W185. doi: 10.1093/nar/gkm321
Morschhauser, J. (2010). Regulation of multidrug resistance in pathogenic fungi. Fungal Genet. Biol. 47, 94–106. doi: 10.1016/j.fgb.2009.08.002
Mortazavi, A., Williams, B. A., Mccue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods. 5, 621–628. doi: 10.1038/nmeth.1226
Murat, C., Payen, T., Noel, B., Kuo, A., Morin, E., Chen, J., et al. (2018). Pezizomycetes genomes reveal the molecular basis of ectomycorrhizal truffle lifestyle. Nat. Ecol. Evol. 2, 1956–1965. doi: 10.1038/s41559-018-0710-4
Park, J., Park, J., Jang, S., Kim, S., Kong, S., Choi, J., et al. (2008). FTFD: an informatics pipeline supporting phylogenomic analysis of fungal transcription factors. Bioinformatics 24, 1024–1025. doi: 10.1093/bioinformatics/btn058
Parra, G., Bradnam, K., and Korf, I. (2007). CEGMA: a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics 23, 1061–1067. doi: 10.1093/bioinformatics/btm071
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Pion, M., Spangenberg, J. E., Simon, A., Bindschedler, S., Flury, C., Chatelain, A., et al. (2013). Bacterial farming by the fungus Morchella crassipes. Proc. R. Soc. B Biol. Sci. 280:9. doi: 10.1098/rspb.2013.2242
Rawlings, N. D., Barrett, A. J., and Bateman, A. (2010). MEROPS: the peptidase database. Nucleic Acids Res. 38, D227–D233.
Rhind, N., Chen, Z., Yassour, M., Thompson, D. A., Haas, B. J., Habib, N., et al. (2011). Comparative functional genomics of the fission yeasts. Science 332, 930–936. doi: 10.1126/science.1203357
Robinson, M. D., Mccarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Ryan, R. P., Germaine, K., Franks, A., Ryan, D. J., and Dowling, D. N. (2008). Bacterial endophytes: recent developments and applications. FEMS Microbiol. Lett. 278, 1–9. doi: 10.1111/j.1574-6968.2007.00918.x
Saier, M. H. Jr., Tran, C. V., and Barabote, R. D. (2006). TCDB: the Transporter Classification Database for membrane transport protein analyses and information. Nucleic Acids Res. 34, D181–D186. doi: 10.1093/nar/gkj001
Sanderson, M. J. (2002). Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Mol. Biol. Evol. 19, 101–109. doi: 10.1093/oxfordjournals.molbev.a003974
Schlosser, A., Meldorf, M., Stumpe, S., Bakker, E. P., and Epstein, W. (1995). TrkH and its homolog, TrkG, determine the specificity and kinetics of cation transport by the Trk system of Escherichia coli. J. Bacteriol. 177, 1908–1910. doi: 10.1128/jb.177.7.1908-1910.1995
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V., and Zdobnov, E. M. (2015). BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212. doi: 10.1093/bioinformatics/btv351
Song, H. Y., Kim, D. H., and Kim, J. M. (2018). Comparative transcriptome analysis of dikaryotic mycelia and mature fruiting bodies in the edible mushroom Lentinula edodes. Sci. Rep. 8:15. doi: 10.1038/s41598-018-27318-z
Spatafora, J. W., and Bushley, K. E. (2015). Phylogenomics and evolution of secondary metabolism in plant-associated fungi. Curr. Opin. Plant Biol. 26, 37–44. doi: 10.1016/j.pbi.2015.05.030
Splivallo, R., Vahdatzadeh, M., Macia-Vicente, J. G., Molinier, V., Peter, M., and Egli, S., et al. (2019). Orchard Conditions and Fruiting Body Characteristics Drive the Microbiome of the Black Truffle Tuber aestivum. Front Microbiol. 10:1437. doi: 10.3389/fmicb.2019.01437
Stamatakis, A. (2006). RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 22, 2688–2690. doi: 10.1093/bioinformatics/btl446
Stanke, M., Diekhans, M., Baertsch, R., and Haussler, D. (2008). Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics 24, 637–644. doi: 10.1093/bioinformatics/btn013
Strasser, K., Mcdonnell, E., Nyaga, C., Wu, M., Wu, S., Almeida, H., et al. (2015). mycoCLAP, the database for characterized lignocellulose-active proteins of fungal origin: resource and text mining curation support. Database (Oxford) 2015:bav008. doi: 10.1093/database/bav008
Taylor, J. W., and Berbee, M. L. (2006). Dating divergences in the Fungal Tree of Life: review and new analyses. Mycologia 98, 838–849. doi: 10.3852/mycologia.98.6.838
Tietel, Z., and Masaphy, S. (2018). True morels (Morchella)—Nutritional and phytochemical composition, health benefits and flavor: a review. Crit. Rev. Food Sci. Nutr. 58, 1888–1901. doi: 10.1080/10408398.2017.1285269
Traeger, S., Altegoer, F., Freitag, M., Gabaldon, T., Kempken, F., Kumar, A., et al. (2013). The genome and development-dependent transcriptomes of Pyronema confluens: a window into fungal evolution. PLoS Genet. 9:e1003820. doi: 10.1371/journal.pgen.1003820
Trappe, M. J., Trappe, J. M., and Bonito, G. M. (2010). Kalapuya brunnea gen. & sp. nov. and its relationship to the other sequestrate genera in Morchellaceae. Mycologia 102, 1058–1065. doi: 10.3852/09-232
Turgeon, B. G., and Yoder, O. C. (2000). Proposed nomenclature for mating type genes of filamentous ascomycetes. Fungal Genet. Biol. 31, 1–5. doi: 10.1006/fgbi.2000.1227
Ullah, A., Nisar, M., Ali, H., Hazrat, A., Hayat, K., Keerio, A. A., et al. (2019). Drought tolerance improvement in plants: an endophytic bacterial approach. Appl. Microbiol. Biotechnol. 103, 7385–7397. doi: 10.1007/s00253-019-10045-4
Vieira, V., Fernandes, A., Barros, L., Glamoclija, J., Ciric, A., Stojkovic, D., et al. (2016). Wild Morchella conica Pers. from different origins: a comparative study of nutritional and bioactive properties. J. Sci. Food Agric. 96, 90–98. doi: 10.1002/jsfa.7063
Wang, M., Gu, B., Huang, J., Jiang, S., Chen, Y., Yin, Y., et al. (2013). Transcriptome and proteome exploration to provide a resource for the study of Agrocybe aegerita. PLoS One 8:e56686. doi: 10.1371/journal.pone.0056686
Wang, Y., Tang, H., Debarry, J. D., Tan, X., Li, J., Wang, X., et al. (2012). MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 40:e49. doi: 10.1093/nar/gkr1293
Wasser, S. P. (2017). Medicinal mushrooms in human clinical studies. Part I. Anticancer, oncoimmunological, and immunomodulatory activities: a review. Int. J. Med. Mushrooms 19, 279–317. doi: 10.1615/intjmedmushrooms.v19.i4.10
Wei, W., Mccusker, J. H., Hyman, R. W., Jones, T., Ning, Y., Cao, Z., et al. (2007). Genome sequencing and comparative analysis of Saccharomyces cerevisiae strain YJM789. Proc. Natl. Acad. Sci. U.S.A. 104, 12825–12830. doi: 10.1073/pnas.0701291104
Wipf, D., Bedell, J. P., Botton, B., Munch, J. C., and Buscot, F. (1996). Polymorphism in morels: isozyme electrophoretic analysis. Can. J. Microbiol. 42, 819–827. doi: 10.1139/m96-103
Xu, J. R., and Hamer, J. E. (1996). MAP kinase and cAMP signaling regulate infection structure formation and pathogenic growth in the rice blast fungus Magnaporthe grisea. Genes Dev. 10, 2696–2706. doi: 10.1101/gad.10.21.2696
Yamada, M., Sakuraba, S., Shibata, K., Taguchi, G., Inatomi, S., Okazaki, M., et al. (2006). Isolation and analysis of genes specifically expressed during fruiting body development in the basidiomycete Flammulina velutipes by fluorescence differential display. FEMS Microbiol. Lett. 254, 165–172. doi: 10.1111/j.1574-6968.2005.00023.x
Yang, J., Wang, L., Ji, X., Feng, Y., Li, X., Zou, C., et al. (2011). Genomic and proteomic analyses of the fungus Arthrobotrys oligospora provide insights into nematode-trap formation. PLoS Pathog. 7:e1002179. doi: 10.1371/journal.ppat.1002179
Yuan, Z. S., Liu, F., Xie, B. G., and Zhang, G. F. (2018). The growth-promoting effects of endophytic bacteria on Phyllostachys edulis. Arch. Microbiol. 200, 921–927. doi: 10.1007/s00203-018-1500-8
Zhang, Q., Liu, W., Cai, Y., Lan, A. F., and Bian, Y. (2018). Validation of internal control genes for quantitative real-time PCR gene expression analysis in Morchella. Molecules 23:2331. doi: 10.3390/molecules23092331
Zhao, Z., Liu, H., Wang, C., and Xu, J. R. (2013). Comparative analysis of fungal genomes reveals different plant cell wall degrading capacity in fungi. BMC Genomics 14:274. doi: 10.1186/1471-2164-14-274
Keywords: Morchella conica SH, genomic, transcriptomic, fruiting body, endophytic bacterial communities
Citation: Lü BB, Wu GG, Sun Y, Zhang LS, Wu X, Jiang W, Li P, Huang YN, Wang JB, Zhao YC, Liu H, Song LL, Mo Q, Pan AH, Yang Y, Long XQ, Cui WD, Zhang C, Wang X and Tang XM (2021) Comparative Transcriptome and Endophytic Bacterial Community Analysis of Morchella conica SH. Front. Microbiol. 12:682356. doi: 10.3389/fmicb.2021.682356
Received: 18 March 2021; Accepted: 10 June 2021;
Published: 20 July 2021.
Edited by:
Hector Mora Montes, University of Guanajuato, MexicoReviewed by:
Li Yu, Jilin Agricultural University, ChinaLi Xiang, Institute of Chinese Materia Medica (CACMS), China
Copyright © 2021 Lü, Wu, Sun, Zhang, Wu, Jiang, Li, Huang, Wang, Zhao, Liu, Song, Mo, Pan, Yang, Long, Cui, Zhang, Wang and Tang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xue M. Tang, c2Fhc194bXRhbmdAZm94bWFpbC5jb20=
†These authors have contributed equally to this work