Microbial communities associated with the black morel Morchella sextelata cultivated in greenhouses

Morels (Morchella spp.) are iconic edible mushrooms with a long history of human consumption. Some microbial taxa are hypothesized to be important in triggering the formation of morel primordia and development of fruiting bodies, thus, there is interest in the microbial ecology of these fungi. To identify and compare fungal and prokaryotic communities in soils where Morchella sextelata is cultivated in outdoor greenhouses, ITS and 16S rDNA high throughput amplicon sequencing and microbiome analyses were performed. Pedobacter, Pseudomonas, Stenotrophomonas, and Flavobacterium were found to comprise the core microbiome of M. sextelata ascocarps. These bacterial taxa were also abundant in the soil beneath growing fruiting bodies. A total of 29 bacterial taxa were found to be statistically associated to Morchella fruiting bodies. Bacterial community network analysis revealed high modularity with some 16S rDNA operational taxonomic unit clusters living in specialized fungal niches (e.g., pileus, stipe). Other fungi dominating the soil mycobiome beneath morels included Morchella, Phialophora, and Mortierella. This research informs understanding of microbial indicators and potential facilitators of Morchella ecology and fruiting body production.


INTRODUCTION
Morels (Morchella spp.) are an iconic genus of edible mushrooms that are distributed across the Northern hemisphere (O'Donnell et al., 2011). Morels have a long history of use in Europe, and are sought after in North America and Asia. They remain an economically important culinary mushroom today, and are commercially harvested in the springtime when they fruit naturally (Obst & Brown, 2000;Pilz et al., 2007). For example, in western North America, morels have been estimated to contribute $5-10 million to the economy through direct sales (Pilz et al., 2007).
Morchella is a species-diverse genus. Classical taxonomic treatments of Morchella based on morphological characters are complicated by the extreme variation in macro-characters in-depth characterizations of fungal and prokaryotic communities associated with M. sextelata and the soils beneath their fruiting bodies.

METHODS
Sampling microbial communities associated with morel fruiting bodies and soils beneath fruiting bodies Morel fruiting bodies and soils beneath growing morels were sampled from a high-tunnel greenhouse in Caohaizi Village, Xundian County, Kunming City, Yunnan Province, China, where the black morel M. sextelata was being cultivated. The site is situated 1,950 m in elevation. The pileus and stipe from five mature (>10 cm) and five immature (<1 cm) fruiting bodies were sampled by placing a piece of tissue roughly one cm 2 in size into CTAB 4X buffer with a flame sterilized razor. Approximately two cm 3 of soil was also sampled from directly below each morel fruiting body. Soils were dried completely with silica beads and were kept on silica until processing (described below). In total, microbial analyses were performed on 20 samples, 10 M. sextelata ascocarps (five young and five mature), and 10 soils beneath the ascocarps, which were analyzed for both bacterial (16S rDNA) and fungal (ITS rDNA) communities. Bacterial communities were determined for 10 morel ascocarps, including pileus (n = 10, five mature and five immature) and stipe (n = 10, five mature and five immature) tissues.

Molecular methods
DNA was extracted from~0.5 g of dried and homogenized soils with the PowerMag Ò Soil DNA Isolation Kit (Qiagen, Carlsbad, CA, USA) following manufacturer's recommendations. Morel tissues were ground with a sterile micro pestle and then extracted using standard chloroform extraction protocol (Trappe, Trappe & Bonito, 2010). Extracted DNA was amplified using DreamTaq Green DNA Polymerase (ThermoFisher Scientific, Waltham, MA, USA) with the following primer sets: ITS1f-ITS4 for Fungi and 515F-806R for Bacteria and Archaea, following a protocol based upon the use of frameshift primers as described by Chen et al. (2018) and Lundberg et al. (2013). PCR products were stained with ethidium bromide, separated through gel electrophoresis, and imaged under UV light. Amplicon concentrations were normalized with the SequalPrep Normalization Plate Kit (ThermoFisher Scientific, Waltham, MA, USA) and pooled. Amplicons were then concentrated 20:1 with Amicon Ultra 0.5 mL 50K filters (EMD Millipore, Darmstadt, Germany) and purified with Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA, USA). A synthetic mock community with 12 taxa and four negative (no DNA added) controls was included to assess sequencing quality (Palmer et al., 2018). Amplicons were then sequenced on an Illumina MiSeq analyzer using the v3 600 cycles kit (Illumina, San Diego, CA, USA). Sequence reads have been submitted to NCBI SRA archive under the accession number PRJNA510627.

Statistical analyses
The otu_table.biom (McDonald et al., 2012) with embedded taxonomy classifications and metadata.txt files for each marker gene were imported into the R statistical environment for analysis (R Core Team, 2018). Before proceeding with the analysis, data were quality filtered to remove OTUs with less than 10 total sequences (Lindahl et al., 2013;Oliver, Callaham & Jumpponen, 2015). OTUs that appeared in the negative controls (i.e., contaminants) were removed across all samples when ≥10 reads were present in any single control. Observed OTU richness (S) (Simpson, 1949), Shannon's diversity index (Hill, 1973), and Evenness (Kindt & Coe, 2005) were used as a-diversity metrics. The Shannon index (H) was calculated as H = −∑ i R =1 p i lnp i where p i the proportion of individuals belonging to the i species in the dataset, while the OTU evenness (E) was calculated as E ¼ H lnðsÞ where H is the Shannon diversity index and S the observed OTU richness. Diversity indexes were with the "specnumber" and "diversity" functions in R package vegan (Oksanen et al., 2019) and with the function "diversityresult" in the package BiodiversityR (Kindt & Coe, 2005). After assessing for data normality and homogeneity of variances significant differences between mean alpha-diversity measures were found with ANOVA and Tukey's tests. Rarefaction curves were used to assess OTU richness from the results of sampling (Figs. S2 and S3). To avoid biases and data loss in some groups of samples due to inherent variations in alpha-diversity in soils compared to morels, OTUs were normalized using the R package metagenomeSeq before calculating β-diversity (Paulson et al., 2013). Principal coordinate analysis (PCoA) was used to investigate community β-diversity with the function "ordinate" from the phyloseq package (McMurdie & Holmes, 2013). Diversity patterns were then tested for statistical differences across sites in the vegan R package with the PERMANOVA function "adonis" and tested for homogeneity of variances with the function "betadisper." OTUs that showed high and significant correlation with sample groups were identified through the function "multipatt" in the indicspecies package (De Cáceres & Legendre, 2009).
To assess co-occurrences among OTUs a bipartite network was produced for the prokaryotic communities with the "spiec.easi" function in the SpiecEasi R package (Kurtz et al., 2015). To build the network, the following parameters were used: lambda.min. ratio=1e-2, nlambda=50, rep.num=99. The network was constructed using the OTUs present in at least 15 samples to increase the sensitivity of the analysis. After assessing network stability using the "getStability" function in SpiecEasi, general (i.e., modularity, sparsity, transitivity) and individual OTUs (i.e., degree, closeness centrality, betweenness centrality, articulation points) network indexes were calculated. The network was visualized with the Fruchterman-Reingold layout with 10 4 permutations as implemented in the igraph R package (Csardi & Nepusz, 2006). A heatmap showing abundances of prokaryotic OTUs statistically associated with Morchella ascocarps was created using the ComplexHeatmap R package (Gu, Eils & Schlesner, 2016). All statistical analyses and graphs were performed in R version 3.4.4 (R Development Core Team, 2018).

High-throughput sequencing results
After quality filtering, a total of 215,201 reads were analyzed with an average read depth of 21,520 across 10 samples for the ITS marker and 2,237,810 reads with an average read depth of 74,593 reads across 30 samples for 16S rDNA. After removing contaminants, as well as negative and mock samples, a total of 509 OTUs for fungal communities and 5,169 OTUs for prokaryotic communities were obtained. Our synthetic mock community matched the 12 artificial taxa, which were sequenced alongside with the samples. No mock sequences were detected in any other libraries indicating that barcode switching was not an issue in this study.

Microbial richness and evenness in soils beneath Morchella fruiting bodies
Significant differences (p ≤ 0.05) in OTU richness of the prokaryotic community were found between soil, stipe, and pileus samples ( Table 1). The soil compartment showed almost 10-fold higher richness than was present in Morchella pileus or stipe compartments. Similar trends were true for both Evenness (E) and Shannon index (H) diversity measurements. No differences were found when average alpha-community measures were compared between young and mature morel samples. Fungal alpha diversity trended to be slightly higher in the samples of the young Morchella, but this was not statistically significant (Table 1).

Fungal and prokaryotic community β-diversity in Morchella samples
Principal coordinate analysis (PCoA) ordination graphs performed on the 16S rDNA data show that the difference between the soil from pileus and stipe prokaryotic communities explained the variance obtained in the first axis (49.9%), while differences between pileus and stipe samples are evident in the second axis (18.3%) ( Fig. 2A). PCoA ordination graphs performed on ITS soil data show that the variance of the first axis (66.3%) is due to differences between samples collected under mature compared to young Morchella fruiting bodies (Fig. 2B). Variation obtained for the second axis (8.8%) is due to the high heterogeneity (See below) of the samples collected under young M. sextelata fruiting bodies. PERMANOVA analysis of the 16S dataset show that there was a significant effect of the maturity stage of Morchella samples on the prokaryotic communities (Table 2). The PERMANOVA analysis of the ITS dataset show that there was a significant effect of maturity stage of M. sextelata fruiting bodies on the soil fungal communities ( Table 2) that was not due to sample group dispersion (Fig. S5).

Indicator species and intersections between stage and site
Several prokaryotic OTUs were significantly associated with the pileus, stipe, or associated soil (and combination of them) portions of M. sextelata fruiting bodies (Table S1). A heatmap of the OTUs associated with Morchella fruiting bodies (i.e., associated to pileus, stipe, stipe and pileus, soil and pileus, stipe and soil) is provided in Fig. 3. Two OTUs were statistically associated to Morchella's pileus: Corynebacterium sp. and Pseudanabaena sp. Two OTUs were also associated to the Morchella stipe: Granulicatella sp. and an unidentified OTU in Coxiellaceae. All other OTUs reported in the heatmap were associated to two different groups. Among these OTUs, one specific Pedobacter sp.1 was  associated to both pileus and stipe and was more abundant in these two compartments than it was in the soils. Venn diagrams show that soil samples contained a high number of unique prokaryotic OTUs (3,239) compared to pileus (63), and stipe (34) samples (Fig. 4A) in contrast to what  was shared among them (789). Most bacterial OTUs detected in Morchella fruiting bodies were found in young and mature fruiting bodies (4,644), with only 194 and 331 uniquely present in young or mature samples, respectively. In the fungal communities, mature and young morel soils shared 396 OTUs, while 33 and 80 OTUs were only present in mature or young specimen, respectively (Fig. 4C).

Network analysis
The bipartite network (140 vertex, 199 edges, stability = 0.044) that was obtained is a sparse network (Figs. 5A and 5B), having a low number of possible edges (sparsity ≈ 2%).   The network showed low transitivity (≈0.15) which is a measure of the probability that the adjacent vertices of a vertex are connected. The network showed high modularity (≈0.5) which measures the division into subgraphs (i.e., communities or modules) in which vertex (i.e., OTUs) are more interconnected together than with the rest of the network. A total of 45 modules were identified, with the first five modules containing 40% of the total OTUs: Module 1 contained 29 OTUs; Module 2 contained nine OTUs; Module 3 contained eight OTUs; Module 4 contained six OTUs; Module 5 contained five OTUs. A total of 17 modules were composed of one single OTU (Fig. 5A). Several modules were peripheral and negatively connected (edge weight max = 0.34, min = −0.20) to other modules. Module 5 contained two indicator OTUs identifies for the pileus and stipe niche. Most of the indicator taxa for the stipe and soil environments were in single OTU modules (see Fig. 3 for taxonomic position), disconnected from the main network. Taxonomic classifications at phylum level of each OTU in the network is shown in Fig. 5B. Proteobacteria, Acidobacteria, and Gemmatimonadetes were dominant in the first five modules (Fig. 5C). Interestingly, archaeal OTUs in the Euryarchaeota, Crenarchaeota were also present in the network. In addition to identifying nodes with high degree (number of connections), some OTUs were identified as articulation points, node whose removal disconnects the network (e.g., OTU_2352).

DISCUSSION
Black morels are cultivated in greenhouse conditions in non-sterilized soils (Liu et al., 2018a). It has been hypothesized that fungi and bacteria living in these substrates may facilitate, or conversely, inhibit developmental transitions and fruiting body development (Liu et al., 2017). Soils where morels are cultivated successfully were highly colonized by Morchella mycelium, especially in soils beneath mature morel fruiting bodies. The morel mycelium inoculated in soils appears to overgrow and potentially exclude other fungal taxa.
Regarding prokaryotic communities, Pedobacter, Pseudomonas, Stenotrophomonas, and Flavobacterium were dominant in the microbiome of M. sextelata fruiting bodies. The high abundance of Pseudomonas (Proteobacteria) in morel fruiting bodies raises questions concerning their roles in the development of morels, following observations on the occurrence and diversity of bacterial communities on Tuber magnatum during truffle maturation, Pseudomonas putida farming by M. crassipes (Pion et al., 2013), and the abundance of Pseudomonas OTUs in soils where black morels are cultivated in Sichuan, China (Liu et al., 2017).
Moreover, the relative abundances of bacterial groups varied between vegetative (stipe) and fertile (pileus) tissues of morel mushrooms, as well as from the soil beneath them. For instance, the pileus of Morchella was enriched in Pseudomonas, Stenotrophomonas, and Flavobacteria compared to stipe microbial communities. The stipe was mostly colonized by Pedobacter (83%) compared to the pileus (39%) and the soil where it accounted for only 0.4% of relative abundance of bacteria. OTUs classified as Pedobacter were statistically associated to pileus and stipe tissues and were present in different modules in the microbial network. This indicates that the pileus tissue may recruit a specific set of prokaryotic taxa which are not recruited to the stipe. This is supported by a significant reduction in prokaryotic richness in the pileus and stipe compared to the soils. Of interest, the two tissue types also smelled different. Previous studies have indicated differences in the chemical composition of Amanita pileus and stipes due to metabolite production in the fruiting body (Deja et al., 2014). If similar chemical differences exist between Morchella pileus and stipe, this could offer an explanation for the existence of different prokaryotic communities within distinct tissues of the Morchella fruiting body and the soil beneath them.
Morchella pileus, stipes, and soils were also shown to be specific niches for other indicator bacterial taxa. Surprisingly, human and animal (sometimes plant) pathogens such as Corynebacterium, Granulicatella, Streptococcus, and Staphylococcus were found exclusively associated to the pileus and/or stipe environment (Collins et al., 2004;Cargill et al., 2012). These taxa are components of the microbial network associated with Morchella fruiting bodies (Fig. 5), although they were found in peripheral modules that were negatively connected with the main structure. Some other taxa such as Lacibacter (Qu et al., 2009) or Sediminibacterium (Qu & Yuan, 2008), which are bacteria common in soil, were also identified as indicator species but were not present in our network.
It has been hypothesized that microbes in the soil are necessary for morel fruiting to occur. The role of Pseudomonas in the cultivation of button mushrooms (Agaricus bisporus) has been studied previously, and was shown to increase both yield and primordia formation (Zarenejad, Yakhchali & Rasooli, 2012;Chen et al., 2013;Pent, Põldmaa & Bahram, 2017). The relative abundance of Pseudomonas species increased throughout cultivation cycle of Agaricus bisporus and peaked around the time of fruiting (Chen et al., 2013). It was also shown that the presence of specific strains of Pseudomonas putida in Agaricus inoculum increased mushroom yields by as much as 14% (Zarenejad, Yakhchali & Rasooli, 2012). Previous research found that Pseudomonas putida stimulates sclerotia formation in Morchella (Pion et al., 2013). These results are consistent with our findings that Pseudomonas are abundant in soils and fruiting bodies of cultivated morels, thus, they may be important in the growth and fruiting of these fungi. Liu et al. (2017) Oh et al. (2018) also demonstrated that Pseudomonas are the most common bacteria overall in soils where morels are cultivated, with the highest abundance in the treatment having the highest yield of morel ascocarps, however, bacterial associated with morel fruiting bodies was not assessed. The effect of Flavobacterium spp. on mushroom fruiting body formation is not well studied, but these bacteria have been shown to be associated with the successful cultivation of Pleurotus ostreatus (Cho et al., 2008). Thus, it is possible that Flavobacterium contribute to the formation of mushroom fruiting bodies. The recruitment of prokaryotic communities by Morchella may occur due to a selection by the fungus for specific taxa, or because it offers a preferential niche for bacterial growth. It is also possible that these two factors act simultaneously. For example, Cantharellus cibarius is populated by millions of different bacteria that are thought to be existing on fungal exudates including trehalose and mannitol (Rangel-Castro, Danell & Pfeffer, 2002). Fast growing bacteria that live on fungal-derived nutrients may occupy this niche quickly and may play a role in inhibiting the entry of other bacteria or pathogens (Liu et al., 2018b). Future studies can directly test these hypotheses by assessing the importance of management and specific bacterial taxa on the morel microbiome and fruiting body production.

CONCLUSIONS
In conclusion, our work adds further evidence that the fungal host plays a role in the selective recruitment of specific bacterial taxa. Our study found that the Morchella microbiome is consistently comprised of a small community of bacteria, including Pedobacter, Pseudomonas, Stenotrophomonas, and Flavobacteria, which appear to be recruited from the soil and enriched in fungal fruiting body tissues. Among those, Pedobacter was enriched in and significantly associated with the pileus environment in respect to the stipe and soil compartments. Although some of the bacteria groups detected on morels have also been detected in other mushrooms, based on this preliminary study, many microbial taxa may be exclusive to Morchella. The role of host identity may provide predictive explanation for differences between microbiomes of morels and other mushrooms. Future research is warranted to test the function of these bacteria on morel fruitification and management.
Program project number 4002 (2018). Reid Longley graduate research fellowship support from the Plant Biotechnology for Health and Sustainability Training Program Project NIH T32-GM110523. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: Michigan State University AgBioResearch NIFA: MICL02416 and GREEEN GR17-083. Science and Technology Service Network Initiative, Chinese Academy of Sciences (2017). Guizhou Science and Technology Program: 4002 (2018). Plant Biotechnology for Health and Sustainability Training Program Project NIH T32-GM110523.