The whole-genome sequence analysis of Morchella sextelata

Morchella are macrofungi and are also called morels, as they exhibit a morel-like upper cap structure. Morels contain abundant essential amino acids, vitamins and biologically active compounds, which provide substantial health benefits. Approximately 80 species of Morchella have been reported, and even more species have been isolated. However, the lack of wild Morchella resources and the difficulties associated with culturing Morchella have caused a shortage in the morels available for daily consumption. Additionally, in-depth genomic and morphological studies are still needed. In this study, to provide genomic data for further investigations of culturing techniques and the biological functions of Morchella sextelata (M. sextelata), de novo genome sequencing was carried out on the Illumina HiSeq. 4000 platform using both the Illumina 150 and PacBio systems. The final estimated genome size of M. sextelata was 52.93 Mb, containing 59 contigs and a GC content of 47.37%. A total of 9,550 protein-coding genes were annotated. In addition, the repeat sequences, gene components and gene functions were analyzed using various databases. Furthermore, the secondary metabolite gene clusters and the predicted structures of their products were analyzed. Finally, a genomic comparison of different species of Morchella was performed.

Second-generation genome sequencing and assembly. Illumina PE150 sequencing generated 2,985 Mb raw reads, and 2,600 Mb clean reads were retained after the quality-control steps (data not shown). Then, the M. sextelata genome size was estimated to be 58.77 Mb after revision by K-mer spectrum analysis with K-mer = 15 and depth 27.36. The heterozygous rate was 0.10%, and the repeat rate was 36.09% (Table S1).
The clean reads were de novo assembled, and their scaffolds and contigs were analyzed. Notably, the generated results were only used to provide a reference for the following PacBio RSII third-generation analysis. Therefore, not all data are shown. third-generation genome sequencing and assembly. The raw reads from PacBio RSII were filtered by removing the reads with low quantity, and 869,182 clean reads were ultimately obtained (Table S2). The mean concordance of reads was 0.84 ( Fig. 1), the mean read length was 6,626 bp and the N50 read length was 10,165 bp (Fig. 2).  Genome component predictions. The interspersed repetitive sequences (IRSs) and the tandem repeats (TRs) were determined. IRSs contain short interspersed nuclear elements (SINEs) and long interspersed nuclear elements (LINEs), with the latter often showing transposition activity. In M. sextelata, a total of 15,485 TRs with repeat sizes of 1-1,894 bp, 10,026 minisatellite DNAs with repeat sizes of 10-60 bp and 2,702 microsatellite DNAs with repeat sizes of 2-6 bp were predicted (Table S3).
Noncoding RNAs (ncRNAs) play essential roles in biological processes, although they do not transcribe proteins. The results showed that there were 445 tRNAs with an average length of 85 bp in the M. sextelata genome.  Table 2. The summary of the final genome size and protein-coding gene annotation.
Gene function. The functions of genes were predicted by different databases. The number of genes enriched in each functional classification is shown in Table 3. Generally, many more coding genes were enriched in Kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) compared to other databases.
In the GO analysis, a total of 5,876 genes were enriched. Among them, 2,184 genes were associated with cell, the same number of genes contributed to cell part, and 1,024 genes were enriched in organelle. A total of 3,381 genes were predicted to have binding abilities. The catalytic activity and transporter activity for 2,986 and 389 genes were predicted. Various biological processes were predicted to be mediated by the genes of M. sextelata. For instance, 3,335, 3,239, and 859 genes possibly participate in cellular process, metabolic process, and localization, respectively (Fig. S3).
KEGG pathway analysis was conducted, and the enrichment of genes in different processing, metabolism, and organismal systems was annotated as shown in Fig. S4. The results showed that 203 genes were enriched in the pathway of infectious disease, 110 and 91 genes were enriched in the signaling of cancers and neurodegenerative disease, and 56 genes were enriched in the process correlated with endocrine and metabolic disease.
The genes were also classified based on COG functional database. The majority of annotated genes (331) were predicted to have a general function only. The top three enriched categories were 1. posttranslational modification, protein turnover, chaperones (240 genes); 2. translation, ribosomal structure and biogenesis (219 genes); and 3. amino acid transport and metabolism (161 genes) (Fig. S5).
NR is a database established by NCBI for the annotation of a wide range of information, including species details. Here, most of the genes were annotated to be nonredundant proteins in Tuber melanosporum, 146 of the genes were annotated in Podospora anserine, and 84 genes were annotated in Pseudogymnoascus pannorum. The top 20 predicted nonredundant protein species are shown in Fig. S6.
Based on the above gene annotation and function analysis, the genomic map was drawn and is shown in Fig. 4. In addition, the function of proteins was classified using TCDB. As shown in Fig. S7, most transporter proteins were enriched in primary active transporters, up to 120. Next, 111 proteins were annotated as electrochemical potential-driven transporters, of which 55 were enriched in channels and pores.
There were 51 genes annotated in the endoplasmic reticulum retrieval protein1 (Rer1; putative heavy metal transporter) family, 30 genes in the major facilitator superfamily (MFS), and 28 genes in the H + -or Na + -translocating NADH dehydrogenase (NDH) family. Additionally, some genes were predicted to be members of transporter families, such as the genes A3238 and A0956 in the ammonia transporter channel (Amt) family, A2965 and A3464 in the CorA metal ion transporter (MIT) family, and A8387, A8696, A9300, A1158, A1940, A4798 in the drug/metabolite transporter (DMT) Superfamily.
The CAZy database includes glycoside hydrolases (GHs), glycosyltransferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), and auxiliary activities (AAs). When the genes of M. sextelata were annotated by the CAZy database, GHs were the most abundant enzymes, as 159 genes were predicted to be GHs. As shown in Table 4, 57 of the genes were predicted to code carbohydrate modules (CBMs), while 44 of the genes were estimated to be AAs. In addition, the functions of M. sextelata genomes were also annotated by the Pfam and Swiss-Prot databases (data not shown).
The secondary metabolite gene clusters were analyzed. A total of 10 clusters, including 4 terpenes, 1 nps, 1 t1pks and 4 other clusters, were identified. The gene numbers for each cluster are marked on the blue bar (Fig. 5, Table 5 and Table 6). Additionally, the products of the gene clusters t1pks (polyketide synthase) and nrps (nonribosomal peptide synthase) were analyzed, and the chemical structures were predicted as shown in Fig. 6.
In the effector analysis, the cytochrome P450 database was BLAST searched to annotate the genes encoding fungal cytochromes. The results showed that 62 of the genes were predicted to be P450s distributed in 6 classes. Specific gene IDs are provided in Table 7.
Then, the shared and specific genes in different species of Morchella were analyzed as shown in Fig. 7. The results showed that M. eximia exhibited the highest specific gene families with 7,347, and M. sextelata also showed 1,139 specific gene families, which is much higher than the number of specific gene families in other species, ranging from 371 to 837. On the other hand, the number of shared gene families was 3,717.
Furthermore, the dispensable genes, which are all genes except the core genes, in each strain were analyzed, and a heatmap is shown in Fig. 8. Consistent with the results in the Venn diagram, M. eximia (QMFK00000000.1) showed a significant difference compared to other Morchella. Additionally, a large number of specific genes were found in M. sextelata. For example, the gene A8276 specifically expressed in M. sextelata was predicted to be involved in bladder cancer, thyroid cancer, and acute myeloid leukemia. Interestingly, the gene A5703 was also functionally related to bladder cancer, thyroid cancer, and acute myeloid leukemia in these Morchella species.

Discussion
The genome sequencing in this study provides the whole-genome sequence of M. sextelata and its related gene annotation for the first time. This study may provide important data for evaluating Morchella species, improving culture techniques by mediating mating behavior, and discovering biologically active compounds. It is important not only for meeting the increasing demand for M. sextelata but also for further research on M. sextelata.
The comparative genomic analysis of M. sextelata and other Morchella species, M. conica, M. importuna, M. septimelata, and M. eximia, was carried out in this study. To provide additional useful information, the gene annotation file generated in this study was also uploaded and may provide useful data for further research on the differences between various Morchella species and their biological functions in the future.
The biological roles of Morchella have been reported in recent years. The potential neuritogenic activity of the aqueous extract from M. importuna 11 , immunomodulatory activity of the exopolysaccharides extracted from M. conica 12 , and antitumor activity of polysaccharides isolated from M. esculenta were determined 4 . In addition, a polysaccharide from M. conica was found to inhibit H 2 O 2 -induced oxidative stress in human embryonic kidney 293 T cells 13 . In addition to these, many other biological functions of the compounds obtained from different Morchella have been identified, suggesting the great value of further studies of Morchella.

Materials and Methods
M. sextelata sample. The self-mated M. sextelata spores were cultured in sandy soil in a plastic greenhouse with controlled light, temperature and humidity in Shagedu town located at the eastern end of the Ordos Plateau in Inner Mongolia, China. The M. sextelata isolate we cultured usually had darker caps and larger body sizes. The average size was approximately 15 cm high. It is typically planted from October to November and is picked up from April to May in the next year. Since Morchella grows at low temperatures, the weather of Inner Mongolia is very suitable.  www.nature.com/scientificreports www.nature.com/scientificreports/ DNA extraction and quality test.
(1) A fresh Morchella sample was collected and ground completely in a mortar and pestle precooled with liquid nitrogen. (2) The sample was then transferred into a 50-ml tube containing 20 ml GP1 and 500 μl mercaptoethanol and mixed thoroughly by shaking in a 55 °C water bath. The sample was centrifuged at 12,500 rpm for 8 min (GP1 buffer: RK116-02). (3) The supernatant was collected, and saturated phenol/chloroform/isoamyl alcohol (25:24:1) was added.
Then, it was centrifuged after mixing at 12,500 rpm for 8 min. (4) The supernatant was collected and mixed with chloroform/isoamyl alcohol. Then, the mixture was centrifuged at 12,500 rpm for 8 min after vortexing. (5) The supernatant was collected again in a 50-ml new tube, and isopropanol was added and mixed well. The sample was placed at −20 °C to allow precipitation for a few minutes. Then, the mixture was centrifuged (12,500 rpm, 8 min) and the supernatant was removed. (6) The pellet was washed twice with 5 ml 75% alcohol. The supernatant was removed carefully, and the pellet was dried. (7) A total of 100 μl ddH 2 O was added to dissolve the DNA sample. Then, 2 μl RNase was added to digest RNA, and the sample was then shaken several times and incubated at 37 °C for 15 min. To ensure that the harvested DNA sample was suitable for further analysis and that there was no contamination, the quality of the harvested DNA was examined by agarose gel electrophoresis.
Genome sequencing and assembly via a second-generation system. First, the DNA was randomly cut into 350 bp fragments by Covaris ultrasonic processor. The fragments were modified with end-repairing and the addition of A tail and sequencing adaptors and were purified and PCR amplified to generate a 350 bp long fragment library. Next, the library was initially quantified by Qubit 2.0, and the insert size of fragments was tested by Agilent 2100. Then, the effective concentration of the library was examined by q-PCR to ensure that the quality of the 350 bp library was fit for further sequencing analysis. Subsequently, the genome of the Morchella sample in the 350 bp libraries was sequenced by the second-generation Illumina PE150 system. The generated reads were filtered by several steps to obtain the clean reads, which were then analyzed through K-mer (K-mer = 15) to estimate the genome size.
The genome for the M. sextelata sample was obtained by assembling the clean data from the Illumina PE150 system by SOAP de novo (version 2.04) 14 , SPAde 15 and ABySS 16 software and subsequently integrating the genome sequence with CISA 17 software. The initially assembled genome was optimized and supplemented by the gapclose (version 1.12) program. Consequently, the result of the assembly was statistically analyzed by scaffolds and contigs after excluding fragments <500 bp. Finally, the GC content and coverage depth were analyzed to determine whether there was sequence contamination and repeat sequences.
All of the results obtained from the Illumina PE150 system were used as the survey for the PacBio RSII third-generation sequencing and further analyses as a source of genomic information.
Genome sequencing and assembly via a third-generation system. After the survey was performed using the Illumina PE150 system, another part of the DNA sample was cut into 20 Kb fragments by Covaris g-TUBE, and the DNA damage and ends of the DNA were repaired. Second, both ends of the fragments were joined with a barcode. Next, these fragments were purified with AM pure PB to generate 20 Kb SMRT Bell libraries. The fragments 20 Kb in length were then screened through BluePippin and purified with AMPure PB beads. Subsequently, the library was quantified by Qubit 2.0, and the insert size of fragments was examined by Agilent 2100. After that, the genome of M. sextelata in the 20 Kb SMRT Bell libraries was sequenced by the PacBio RSII third-generation sequencing system. The obtained raw reads were filtered to obtain the clean reads.
The clean reads were assembled using SMRT Link v5.1.0 software (https://www.pacb.com/support/ software-downloads/) 18,19 by consulting the results from the second-generation sequencing data. De novo assembly was carried out based on the results from the Illumina PE150 system since no reference genome of other Morchella was available. Then, the genome obtained from the PacBio RSII third-generation sequencing system was used for all of the further analyses, including the prediction of genome size, genome annotation, genome component and gene function.
Such an assembly that combines second-and third-generation sequencing data can revise and decrease the interference of abnormal GC contents, high repetition and hybridity. Therefore, it can improve the integrity and uniformity of the generated genome sequence. It is particularly applicable for complex fungi such as Morchella, which lack comparable genomic information, such as reference sequences, GC contents, and repetitive sequences.
Both sequencing procedures were performed using the Illumina HiSeq. 4000 sequencing platform at Beijing Novogene Bioinformatics Technology Co., Ltd.
Genome annotations. For M. sextelata, by default, we performed de novo prediction using the Augustus 2.7 program 20 to obtain the protein-coding genes.
Genome component predictions. Genome component prediction included the prediction of the coding genes, the repetitive sequences and the ncRNA. The interspersed repetitive sequences were predicted using