- Open Access
Influence of Agaricus bisporus establishment and fungicidal treatments on casing soil metataxonomy during mushroom cultivation
BMC Genomics volume 23, Article number: 442 (2022)
The cultivation of edible mushroom is an emerging sector with a potential yet to be discovered. Unlike plants, it is a less developed agriculture where many studies are lacking to optimize the cultivation. In this work we have employed high-throughput techniques by next generation sequencing to screen the microbial structure of casing soil employed in mushroom cultivation (Agaricus bisporus) while sequencing V3-V4 of the 16S rRNA gene for bacteria and the ITS2 region of rRNA for. In addition, the microbiota dynamics and evolution (bacterial and fungal communities) in peat-based casing along the process of incubation of A. bisporus have been studied, while comparing the effect of fungicide treatment (chlorothalonil and metrafenone). Statistically significant changes in populations of bacteria and fungi were observed. Microbial composition differed significantly based on incubation day, changing radically from the original communities in the raw material to a specific microbial composition driven by the A. bisporus mycelium growth. Chlorothalonil treatment seems to delay casing colonization by A. bisporus. Proteobacteria and Bacteroidota appeared as the most dominant bacterial phyla. We observed a great change in the structure of the bacteria populations between day 0 and the following days. Fungi populations changed more gradually, with A. bisporus displacing the rest of the species as the cultivation cycle progresses. A better understanding of the microbial communities in the casing will hopefully allow us to increase the biological efficiency of the crop.
Mushrooms have been cultivated since ancient times especially in the eastern countries . They have been appreciated not just for their flavour but for their nutritional and medicinal value. About 7000 species of mushrooms are believed to have varying degrees of edibility, and more than 3000 species of 231 different genera are considered to be major edible mushrooms . Despite this, only about 60 edible mushrooms are grown commercially, and around 10 are produced on an industrial scale. Agaricus bisporus (J. E. Lange) Imbach is cultivated all over the world to produce button mushrooms. These carpophores have relative high protein content and contain fiber, vitamins, minerals and bioactive compounds. A. bisporus is grown in compost consisting of wheat straw, horse or chicken manure . After the compost is colonized, it is covered with a casing soil, which is required for maintaining moisture, developing a specific microbial flora and trapping several metabolites [4,5,6]. This casing soil is prepared with black or blonde peat and calcium carbonate, to correct pH 7–8. Peat is a fossil organic material, coming from deposits of vegetal organic matter (moss) accumulated over thousands of years. The use of peat in A. bisporus cultivation is essential, since until now no other material is known to have its characteristics: high water-holding capacity but water that is available for mushroom use, heavy structure but porous to facilitate gas exchange with a specific microbiota . Peat used for mushroom growing accounts for only a small portion of the total volumes of peat extracted worldwide. However, peat is an unsustainably sourced non-renewable resource whose availability is diminishing due to reserve depletion and regulatory restrictions. Many attempts have been made to replace peat by other materials such as pine bark or coconut fibre. However, due to the particular physico-chemical properties and native microbiota, peat is still seen as the only optimal option . Casing soil characteristics may also affect the development of mushroom diseases, such as dry bubble caused by Lecanicillium fungicola (Preuss), Zare and Gams, and wet bubble caused by Mycogone perniciosa (Magnus) Delacroix [9, 10]. To prevent disease occurrence, chemical products are applied. Notably, the approved fungicides for mushroom use were not specifically developed for mushroom pathogens, and therefore the mushroom industry face side problems such as reduced specificity, resistance outbreaks among pathogens and regulatory limitations for their use. Currently the number of permitted active ingredients with fungicidal action is very limited, therefore it is important to investigate in this line to find alternatives in the fight against these pathogens.
To improve mushroom cultivation methods, the role of microbial communities in the casing soil must be understood. The knowledge generated may contribute to the development of alternative materials of non-fossil origin, more environmentally friendly and with inhibiting properties of other pathogenic organisms, thus improving both the cultivated mushroom sector and other agricultural sectors that use peat, as nurseries or gardens.
In this work we have employed high-throughput techniques by next generation sequencing (NGS) to screen the microbiota dynamics and evolution (bacterial and fungal communities) in peat-based casing soil along the process of fructification of A. bisporus, and the impact of fungicide treatment (chlorothalonil and metrafenone) on the microbial communities. To reach these objectives two highly conserved genetic regions were sequenced in an Illumina MiSeq: V3-V4 of the 16S rRNA gene for bacteria and the ITS2 region of rRNA for fungi.
Data from the trial were divided in three data-set. The first one (data-set DAYS) represents the differences between incubation days in the control without fungicides application: casing day (D00), fungicides application (D04), ruffling day (D07) and beginning of first flush (D16). Ruffling is an agronomic technique used in mushroom farming that consists of raking the surface of the casing soil in order to aerate it and mix compost with the casing so that mushroom colonization is faster and more homogeneous. The second (data-set TREATMENTS D07) represents the differences between the control and the two fungicide treatments (chlorothalonil and metrafenone) at ruffling day (D07). And finally, the third (data-set TREATMENTS D16) represents the differences between the control and the two fungicide treatments in the beginning of the first flush (D16). For these three data-set both bacterial and fungi communities have been analysed.
Sequencing and assembly
Bacterial 16S analysis
For both, the Bacterial 16S analysis and the Fungi ITS2 analysis, all samples (n = 24) were sequenced on the Illumina MiSeq. Sample 26 (G9) from the Bacterial 16S analysis, was excluded from the analysis because of the low number of reads obtained (1168 reads). All data from the metataxonomic analysis are collected in Table 1.
Alpha rarefaction curves were analysed (Supplementary material Figs. 7 and 8), to confirm that sampling depth was enough to observe the full community diversity. Alpha and beta diversity box-and-whisker plots and statistics are collected in Supplementary material Figs. 1, 2, 3, 4, 5 and 6.
Analysis data-set DAYS
First, we performed a metataxomonic analysis to study the different bacterial populations. Phylum and family level were visualized in barplots. Then, we analysed the diversity to check if there were significant differences within samples (alpha diversity) and between samples (beta diversity). Finally, we carried out an ANCOM test to identify different abundances across bacterial samples at genus level.
The most highly represented bacteria phylum was Proteobacteria and Bacteroidota for all samples in data-set DAYS (Fig. 1A). Proteobacteria at D00 appears at 40.34%+/− 1.39 SD (percentage of ASVs classified) and increases along days, with a percentage of 53.70%+/− 0.96 SD at D16. Firminutes phylum presents a 10.88%+/− 0.88 SD the first day (D00), decreasing along the cycle until it reaches a percentage of 1.24%+/− 0.01 SD D16. Other phyla present in smaller percentages, such as Patescibacteria, Verrucomicrobiota, Chloroflexi and Planctomycetota, increased during the days. One candidate phylum (WPS-2 = Eremiobacterota) and one phylum without identification was found.
At family level differences between days can be observed (Fig. 1B). Day D00 shows singularities comparing to the other three days. Most abundant family at day D00 was Bacteroidaceae, with a 6.70%+/− 1.03 SD. This family drastically decreased throughout the crop cycle until D16 with a 0.13% +/− 0.11 SD. Similarly, family Lachnospiraceae almost disappeared at the end of the cycle. Family Xanthomonadaceae was present in all samples with a similar percentage (4–6%). However, the presence of some families increased during days, such as Sphingomonadaceae (from 1.36% +/− 0.47 SD to 7.99% +/− 0.70 SD), Chitinophagaceae (from 1.45% +/− 0.34 SD to 5.03% +/− 0.14 SD), Devosiaceae (from 0.94% +/− 0.07 SD to 2.86% +/− 0.74 SD) or Spirosomaceae (0.20% +/− 0.06 SD to 2.25% +/− 0.20 SD).
Alpha diversity results showed that there were not significant differences among days according to Faith’s phylogenetic diversity (p = 0.059) and Pielou’s species evenness (p = 0.053). Bacterial community richness and evenness showed significant differences by Shannon’s index (p = 0.027). However, results showed higher diversity with-in samples at day 4 (D04) and lower at day zero (D00). Beta diversity results showed significant differences by Bray Curtis: abundance without phylogeny (p = 0.001) and unweighted UniFrac distance: presence and absence of OTUs (Operational Transcriptomic Unit) with phylogeny (p = 0.001).
An ANCOM test was also carried out, to compare the composition of the microbiome within bacterial populations. Results identified 19 different abundances profiles across samples by days at genus level. Aquicella, Halomonas, Solitalea, Lacihabitans, Caulobacter and Anaerosporobacter had the highest W value (Supplementary material Table S1).
The most frequently identified fungal phyla are Ascomycota and Basidiomycota. Noteworthy, significant number of sequences were not classified in any phylum (Fig. 2A). At day zero (D00) the most abundant phylum was Ascomycota (69.13%+/− 0.27 SD). This trend changes along the cultivation cycle with phylum Ascomycota decreasing D04 (57.90%+/− 2.87 SD), D07 (18.68%+/− 7.22 SD) until D16 with 3.86%+/− 1.7 SD. Phylum Basidiomycota increases from D00 (22.38%+/− 0.10 SD) until D16 (91.95% +/− 3.29 SD), becoming the dominant phylum from D07. This trend can be also observed at Family level (Fig. 2B), the Agaricaceae increasing from D00 (2.44%+/− 0.73 SD) to the D16 (89.43%+/− 4.33 SD).
At family level the change in the populations of the different fungi families has been clearly observed. Family Agaricaceae displaced the other families along the crop cycle while A. bisporus mycelium is growing in the casing material. In the first samples, Helotiaceae family was present in an important percentage (18.74%+/− 4.29 SD) that almost disappears at D16 (0.54%+/− 0.13 SD), as well as other families such as Piskurozumaceae, Exidiaceae or Myxotrichaceae. On the other hand, the presence of the family Aspergillaceae increased from D00 (13.54%+/− 6.32 SD) to the D04 (28.84%+/− 1.23 SD) and then decreased by day D07 and D16 (1.58%+/− 0.94 SD).
Alpha diversity results showed significant differences among days using Faith’s diversity (p = 0.024), Pielou’s evenness (p = 0.023) and Shannon’s index (p = 0.024). Beta diversity results showed also significant differences by Bray Curtis: abundance without phylogeny (p = 0.001) and unweighted UniFrac distance (p = 0.001).
ANCOM test results identified different abundance across samples by days at genus level in genus Agaricus (Supplementary material Table S2).
Analysis data-set TREATMENTS D07
Two different fungicide treatments were applied in the trial on day D04: chlorothalonil and metrafenone. Casing samples taken from a control without any treatment were also analysed. The most representative bacterial phyla were Proteobacteria and Bacteroidota for all samples in the data-set TREATMENTS D07 (Fig. 3A). The highest Proteobacteria abundance corresponds to chlorothalonil treatment with 48.51%+/− 1.84 SD. Three candidate phyla were found, WPS-2, WS2 and FCPU426. Most abundant families are common in the treatments and the control (Fig. 3B), such as Flavobacteriaceae (8–9%), Xanthomonadaceae (6.7–6.9%) and Sphingomonadaceae (6.8–5.3%). Family Chitinophagaceae showed higher abundance in metrafenone treatment (5.4% +/−.036 SD) and lower in chlorothalonil treatment (4.51% +/− 0.32 SD). Other families such as Rhizobiaceae had a higher abundance in the control (2.64% +/− 0.43 SD) than with the treatments (chlorothalonil 1.90 +/− 0.70 SD and Metrafenone 1.02% +/− 0.54 SD). Family Rhodanobacteraceae presence is 1% higher with Metrafenone treatment 2.42% +/− 0.13 SD.
Alpha diversity results showed that there were not significant differences between all treatments using Faith’s diversity (p = 0.956), Pielou’s evenness (p = 0.875) and Shannon’s index (p = 0.956). Beta diversity results showed no significant differences by Bray Curtis: abundance without phylogeny (p = 0.149) and unweighted UniFrac distance (p = 0.237).
ANCOM test results did not identify different abundances across samples by treatments.
For the fungal phylum level (Fig. 4A) in the control without treatment, Basidiomycota was the most abundant phylum with 69.71%+/− 10.22 SD. For chlorothalonil treatment both phyla were in similar proportions: Basidiomycota 38.55%+/− 4.95 SD and Ascomycota 38.33%+/− 11.88 SD. For this treatment phylum Mucoromycota was more abundant than in the other 2 treatments. For metrafenone treatment the most abundant phylum was Basidiomycota with 68.43%+/− 7.20 SD.
At family level there were around 10–20% of sequences unidentified (Fig. 4B). Family Agaricaceae in the control and with metrafenone treatment had a high presence (63.48%+/− 11.15 SD and 61.88%+/− 9.41 SD) while in the chlorothalonil treatment the percentage of this family was just 26.09%+/− 8.04 SD. The other two families identified with an important presence were Aspergillaceae and Piskurozymaceae.
Alpha diversity results showed non-significant differences between all treatments using Faith’s diversity (p = 0.066), Pielou’s evenness (p = 0.060) and Shannon’s index (p = 0.060). Beta diversity results showed significant differences by Bray Curtis: abundance without phylogeny (p = 0.013) and unweighted UniFrac distance (p = 0.005).
ANCOM test results do not identify different abundance across samples by treatment at genus level.
Analysis data-set TREATMENTS D16
Similar phylum abundance was found for the data-set TREAMENTS D16 (Fig. 5A). Proteobacteria phylum also showed the highest presence in all samples, representing more than 50% of the bacteria identified. Likewise, for this data-set, Proteobacteria was more abundant in chlorothalonil treatment with 54.78%+/− 2.34 SD. At family level (Fig. 5B) the most common families in all samples were Sphingomonadaceae (8.69%+/− 0.84 SD at control, 9.03%+/− 2.26 SD with chlorothalonil and 4.54%+/− 1.29 SD with metrafenone treatment), Chitinophagaceae, Comamonadaceae, Flavobacteriaceae and Rhizobiaceae.
Alpha diversity showed non- significant differences between all treatments using Faith’s diversity (p = 0.235), Pielou’s evenness (p = 0.367) and Shannon’s index (p = 0.986). Beta diversity results showed significant differences by Bray Curtis: abundance without phylogeny (p = 0.019) and unweighted UniFrac distance (p = 0.011).
ANCOM test results identified different abundance across samples by treatments at genus level. There were 5 genera identified: “uncultured”, Pirellula, Vicingus, 37–13 and Sandaracinus (Supplementary material Table S3).
In this data-set the fungal phylum Basidiomycota was the most abundant in all treatments (Fig. 6A). However, differences in the abundance could be observed. In the case of metrafenone treatment, Basidiomycota represented 70.72%+/− 22.34 SD and Ascomycota at 15.35%+/− 15.29 SD. These standard deviations were higher than in the other treatments. For chlorothalonil treatment, Basidiomycota phylum had a presence of 76.29%+/− 7.91 SD and for the control, Basidiomycota was the most abundant with 91.95%+/− 3.29 SD.
At family level Agaricaceae was the dominant family (Fig. 6B). However, its presence was lower in the treatments (chlorothalonil 69.37%+/− 13.65SD and metrafenone 63.28%+/− 28.38 SD) comparing to the control (89.43%+/− 4.33 SD). In this data-set there were also a high number of unidentified sequences. Similarly to the data-set TREATMENT D07, other families with an important presence were Aspergillaceae and Piskurozymaceae.
Alpha diversity results showed no-significant differences between all treatments using Faith’s diversity (p = 0.670), Pielou’s evenness (p = 0.193) and Shannon’s index (p = 0.288). Beta diversity results showed no-significant differences by Bray Curtis: abundance without phylogeny (p = 0.324) and unweighted UniFrac distance (p = 0.641).
ANCOM test does not identify differences in abundance across samples by treatment at genus level.
The goal of our research is to understand the changes in the microbial dynamics that occur during mycelium development and mushroom (A. bisporus) establishment in the casing soil, and the influence of two fungicide treatment applications on the microbial population changes of this active ecosystem.
The fungicide treatments used in this study were chlorothalonil and metrafenone. Chlorothalonil is a wide spectrum fungicide, which belongs to the chloronitrile substance group. It has been used in mushroom industry mainly against Mycogone perniciosa and Lecanicillium fungicola. Nevertheless, it has been reported to be effective against Cladobotryum mycophilum [11, 12] and Trichoderma spp. , although the selectivity of this active substance is reduced . However, chlorothalonil resistance in Lecanilicium fungicola has been reported . The use of chorothalonil is no longer allowed in Europe, but it is in North America, together with thiabendazole and thiophanate-methyl . In Europe just metrafenone and prochloraz are approved for application in white button mushroom production . Metrafenone belongs to the substance group of benzophenone and is used primarily in cereal production. The mode of action of metrafenone interferes with hyphal morphogenesis , and in mushroom industry has been used against Cladobotryum mycophilum .
The structure of the populations of both bacteria and fungi, at the beginning of the crop cycle (day zero = D00), fluctuate with respect to the following days. As reported in our previous project  data analysed confirm that the casing soil in contact with the compost colonized by A. bisporus, undergoes changes in the dynamics of its population structures. In agreement with other research on A. bisporus, we observed that the predominant bacterial phylum identified in all the casing soil samples was Proteobacteria, followed by Bacteroidota [19, 20]. Proteobacteria population increased along the days as seen in other studies , and Firmicutes phylum decreased throughout the cycle. Other phyla present in smaller percentages, such as Patescibacteria, Verrucomicrobiota, Chloroflexi and Planctomycetota, increased during the days. Most abundant family at day D00 was Bacteroidaceae, which decreased deeply during the cycle, like family Lachnospiraceae that almost disappeared at the end of the cycle. However, the presence of some families increased during days, such as Sphingomonadaceae, Chitinophagaceae, Devosiaceae or Spirosomaceae. Results identified different abundances at genus level across samples by days in 19 genera. The genera with highest W value were Aquicella, Halomonas, Solitalea, Lacihabitans, Caulobacter and Anaerosporobacter. These genera were not present at day D00 and they started to grow to significant populations during the incubation. Between day D07 and D16 a change in family’s abundance was observed. Flavocteriaceae and Xanthomonadaceae had a higher presence at day D07 while Sphingomonadaceae and Chitinophagaceae at day D16. At genus level, Flavobacterium was the most abundant in all samples and Pseudomonas increase along days, as has been already published in our previous research . Five genera with different abundance have been identified between treatments at D16: “uncultured”, Pirellula, Vicingus, 37–13 (unclassified Bacteroidetes) and Sandaracinus. Genus Pirellula seems to be dominant during pinning  and during the first flush in compost samples . These results also confirm the chlorothalonil effect on soil bacterial , on several groups such as Flavobacterium (Vicingus) and δ-proteobacteria (Sandaracinus). The most abundant bacterial species identified (Supplementary material Fig. 9) at day D00 was Bacteroides graminisolvens, a xylanolytic anaerobe bacterial strain  that could be related with peat formation in natural peatland since it takes place in waterlogged soils with little or no access to oxygen. This genus almost disappears along the incubation days.
At day D00 the most abundant fungal phylum was Ascomycota, but throughout the cultivation cycle, the phylum Basidiomycota became dominant as described previoustly [26, 27]. As aforementioned, A. bisporus replaced the other species and colonized the casing soil material. For this reason, from day D16 it becomes very difficult to determine the presence or changes in the populations of other fungi, since reads other than Agaricus/Basidiomycota are residual. Families Helotiaceae, Piskurozymaceae, Exidiaceae or Myxotrichaceae were present in an important percentage at D00 and almost disappeared at D16. Other studies found abundance of genus Agaricus, Apiotrichum, Meliniomyces, Mycothermus, Candida and Pseudeurotium . At species level (Supplementary material Fig. 12) we found the presence of Agaricus, Apiotrichum, Meliniomyces and Candida, among others, but not Meliniomyces or Mycothermus. On the other hand, Aspergillaceae family presence increased from day D00 to day D04 and then it decreased at days D07 and day D16. At day D07 for the control without treatment, Basidiomycota was the most abundant phylum while for chlorothalonil treatment both Basidiomycota and Ascomycota showed similar abundance. In chlorothalonil treatment, phylum Mucoromycota was more abundant than in control and metrafenone. In metrafenone treatment, the most abundant phylum was Basidiomycota. Family Agaricaceae in control and with Metrafenone treatment had a high presence while in chlorothalonil treatment the presence of this family is lower, something that could be associated to the reduced selectivity of this active substance which provokes a delay in the colonisation of the casing by the host . The other two families identified with an important presence were Aspergillaceae and Piskurozymaceae. At D07 we observed an effect of chlorothalonil treatment on the development of A. bisporus, being lower its presence in the chlorothalonil treated samples than in the others (Supplementary material Fig. 13). The toxic effect of chlorothalonil to mushrooms when incorporated into the casing has been reported since 1978 and is confirmed by our results [12, 22, 24]. Beta diversity results corroborated fungicide treatments effect on the bacterial and fungi populations. At D16 phylum Basidiomycota was the most abundant in all treatments. Metrafenone treatment had higher standard deviations than other treatments. At family level Agaricaceae was the dominant family, with a remarkable higher presence in the control. Reads associated to the presence of fungal diseases were detected in the samples as reported in previous works [18, 26]. Lecanicillium genus was detected in samples at D00 and D04, but also in one sample at D07, treated with chlorothalonil. A small number of sequences belonging to Mycogone perniciosa were found just in one sample at D00. The presence of mycoparasites in raw materials (D00) can be associated to a contamination prior to the application in crop. This is very relevant since an early infection can provoke important damages due to heavy disease symptoms.
For the metataxonomy analysis with fungal ITS2, many sequences were assigned to “unidentified” taxon. These lack of results was also found in other similar studies [22, 29] although they were focused on the compost dynamics and not the casing soil changes. In the casing layer unidentified OTUs were also detected in our previous study, mostly in the raw material . Unfortunately, through metataxonomic analysis is not always possible to reach species level, only some sequences of bacteria and fungi have been assigned to the level of species. There is still much work to do in taxonomic identification to build better and more complete databases.
The aim of the study was to understand the changes experienced by different populations of fungi and bacteria in the casing soil during button mushroom production in order to improve peat microbial knowledge and find more sustainable alternatives to peat-based casing. This study shows that the use of QIIME2 is valid for NGS sequence analysis in this type of samples. Noteworthy, available resources in databases do not allow metataxonomic studies to reach the species level in most cases, just the genus level which is a limitation for the depth of the analysis. Significant differences have been found between the populations of both fungi and bacteria between the days of mycelium incubation, being more pronounced between the casing day (D00) and day (D04). Differences have also been found between the fungicide treatments applied. Understanding the effect of chemical treatments on the microbiome equilibrium in the casing can help to design more specific and safe fungicides to cope with mushroom pathogens. More detailed studies are also required to explore the relationship between microbial activity and diversity in casing soil and its implications in the mycelium and primordia development. Finally, a more thorough understanding of the pathogens control would have the potential to increase the quality and quantity of mushrooms produced.
Cultivation room was filled with phase III compost from a commercial mushroom compost yard (Germinados de Lodosa, L.S., Lodosa, Spain). The A. bisporus strain used was Sylvan A15. The casing soil used was 50% black peat and 50% blond peat from Valimex S.L (Andosilla, Spain). Casing pH was adjusted to 8 with limestone addition by the provider. Three biological replicates of the casing samples were destructively taken from the cultivation room along the casing incubation process, from the peat-based casing at four different days (Table 2): casing day (D00), fungicides application day (D04), ruffling day (D07) and at beginning of first flush (D16). A randomized crop design in which fungicide treatments were applied in separated 1m2 areas was set up (three replicates per treatment). Chlorothalonil (Banko Champiñón, chlorothalonil 50%, Arysta Lifescience Spain) treatment concentration was 2 ml/m2 and metrafenone (Vivando, metrafenone 50%, Basf) treatment application was 0.5 ml/m2.
Total DNA extraction
Three biological replicates of the casing samples were studied by extracting genomic DNA (n = 3 replicates per sample type). Fresh samples were homogenized in a ceramic mortar with liquid nitrogen. DNA was extracted from up to 0.5 g of casing with NucleoSpin® Soil kit (MACHEREY-NAGEL). DNA quantity and quality were checked using 2 μl of the purified template in a Qubit 2.0 Fluorometer and the Qubit dsDNA BR Assay kit (Thermo Fisher Scientific, MA, USA) and finally it was visualized on a 1.5% agarose gel stained with Midori-Green Advance staining (Nippon Genetics, Tokyo, Japan).
Library preparation for 16S rRNA and ITS2 gene amplicon sequencing was performed separately following the Illumina (San Diego, CA, USA) recommendations with some modifications. A 2-step amplification procedure was used  using paired end universal bacterial primers  for the V3-V4 hypervariable region of the 16S rRNA gene. For the hypervariable region ITS2, the primers ITS3 and ITS4 were used . The primers used contained the Illumina sequencing adapters (overhang nucleotide sequences) added to the gene-specific sequences as described in previous research .
Sequencing data were demultiplexed using Illumina bcl2fastq© program. Demultiplexed paired FASTQ sequences were imported in QIIME2 artifact format and analysed with QIIME2 v2020.8. Used workflow is described in detail in the following link: https://github.com/Marylou8/Metataxonic-analysis-using-Qiime2-workflow. Quality control was carried out using the DADA2 pipeline  incorporated into QIIME2 [34, 35]. The DADA2 program filtered out PhiX reads, removed chimeric sequences and assigned reads into Amplicon Sequence Variants (ASVs). Taxonomic annotation for bacteria was obtained using SILVA v138 database . Taxonomic annotation for fungi was obtained using UNITE v8.2 2020 database . Chloroplast and mitochondrial contaminants were detected and filtered using the QIIME2 “taxa filter-table” and “taxa filter-seqs” commands.
To filter out low-abundance features, we follow the approach of Morton (https://forum.qiime2.org/t/ancom-giving-strange-w-values/1002/11) where those features which do not sum at least 10 sequences among all samples, as well as those that only appear in one sample were filter out (command “feature-table filter-features”). Differential abundance analysis was analysed with ANCOM test . The data-set will be submitted to National Center for Biotechnology Information (NCBI).
Sequencing statistical analyses were done using QIIME2 v2020.8 . Alpha diversity (within-samples) was analysed using Faith’s phylogenetic diversity: bacterial community richness that incorporates phylogenetic relationships between taxa , Pielou’s species evenness: bacterial community evenness  and Shannon’s index: bacterial community richness and evenness . To determine significance in alpha diversity, non-parametric Kruskal-Wallis comparisons were performed . Box-and-whisker plots for species richness and evenness were generated using QIIME2. Alpha rarefaction curves were analysed, to assess if sampling depth was enough to observe the full community diversity . Beta diversity was analysed using Bray–Curtis distance (abundance without phylogeny)  and unweighted UniFrac distance (presence and absence of OTUs with phylogeny) . Principle Coordinate Analysis (PCoA) plots were generated from Bray-Curtis distances and unweighted UniFrac distance using QIIME2. To generate taxonomy heatmaps R-Studio Version 1.3.1093 was used . Two libraries were necessary to plot the heatmaps from QIIME2 artifacts “tidyverse”  and “qiime2R” .
Availability of data and materials
The datasets analysed during the current study are available in the NCBI SRA repository, with the BioProject accession number PRJNA772891, https://www.ncbi.nlm.nih.gov/sra/PRJNA772891.
Amplicon Sequence Variant
Analysis of Composition of Microbiomes
Divisive Amplicon Denoising Algorithm
National Center for Biotechnology Information
Next Generation Sequencing
Operational Transcriptomic Unit
Polymerase chain reaction
Principal Coordinates Analysis
Quantitative Insights Into Microbial Ecology version 2
Chang ST, Hayes WA. The biology and cultivation of edible mushrooms. New York: Academic; 1976.
Chang ST, Wasser SP. The role of culinary-medicinal mushrooms on human welfare with a pyramid model for human health. Int J Med Mushrooms. 2012;14(2):95–134.
Pelkmans JF, Vos AM, Scholtmeijer K, Hendrix E, Baars JJP, Gehrmann T, et al. The transcriptional regulator c2h2 accelerates mushroom formation in Agaricus bisporus. Appl Microbiol Biotechnol. 2016;100(16):7151–9.
Bels-Koning H. Experiments with casing soils, water supply and climate. Mushroom Sci. 1950;1:78–84.
Flegg P. The casing layer in the cultivation of the mushroom (Psalliota hortensis). J Soil Sci. 1956;7:168–76.
Kalberer PP. Water potentials of casing and substrate and osmotic potentials of fruit bodies of Agaricus bisporus. Sci Hortic. 1987;32(3-4):175–82.
Pardo-Giménez A, Pardo JE, Zied DC. Casing materials and techniques in Agaricus bisporus cultivation. In: Diego CZ, Pardo-Giménez A (eds) Edible and medicinal mushrooms: technology and applications. Chichester: Wiley; 2017.
Carrasco J, Tello ML, Perez M, Preston G. Biotechnological requirements for the commercial cultivation of macrofungi: substrate and casing layer. Chapter 7. In: Singh BP, Chhakchhuak L, editors. Biology of macrofungi. Berlin: Springer Nature; 2018. p. 159–75.
Berendsen RL, Kalkhove SIC, Lugones LG, Wösten HAB, Bakker PAHM. Germination of Lecanicillium fungicola in the mycosphere of Agaricus bisporus. Environ Microbiol Rep. 2012;4(2):227–33.
Largeteau ML, Savoie J-M. Microbially induced diseases of Agaricus bisporus: biochemical mechanisms and impact on commercial mushroom production. Appl Microbiol Biotechnol. 2010;86(1):63–73.
Potočnik I, Vukojević J, Stajić M, Rekanović E, Milijašević S, Stepanović M, et al. Toxicity of fungicides with different modes of action to Cladobotryum dendroides and Agaricus bisporus. J Environ Sci Heal - Part B Pestic Food Contam Agric Wastes. 2009;44(8):823–7.
Carrasco J, Navarro MJ, Santos M, Gea FJ. Effect of five fungicides with different modes of action on cobweb disease (Cladobotryum mycophilum) and mushroom yield. Ann Appl Biol. 2017;171(1):62–9.
Kosanović D, Potočnik I, Vukojević J, Stajić M, Rekanović E, Stepanović M, et al. Fungicide sensitivity of Trichoderma spp. from Agaricus bisporus farms in Serbia. J Environ Sci Heal - Part B Pestic Food Contam Agric Wastes. 2015;50(8):607–13.
Bonnen AM, Hopkins C. Fungicide resistance and population variation in Verticillium fungicola, a pathogen of the button mushroom, Agaricus bisporus. Mycol Res. 1997;101(1):89–96.
Stanojević O, Berić T, Potočnik I, Rekanović E, Stanković S, Milijašević-Marčić S. Biological control of green mould and dry bubble diseases of cultivated mushroom (Agaricus bisporus L.) by Bacillus spp. Crop Prot. 2019;126(September):104944.
Grogan HM, Gaze RH. Fungicide resistance among Cladobotryum spp. - causal agents of cobweb disease of the edible mushroom Agaricus bisporus. Mycol Res. 2000;104(3):357–64.
Opalski KS, Tresch S, Kogel KH, Grossmann K, Köhle K. Hückelhoven R Metrafenone: studies on the mode of action of a novel cereal powdery mildew fungicide. Pest Manag Sci. 2006;62:393–401.
Carrasco J, Tello ML, De Toro M, Tkacz A, Poole P, Péerez-Clavijo M, et al. Casing microbiome dynamics during button mushroom cultivation: implications for dry and wet bubble diseases. Microbiol (United Kingdom). 2019;165(6):611–24.
Siyoum NA, Surridge K, van der Linde EJ, Korsten L. Microbial succession in white button mushroom production systems from compost and casing to a marketable packed product. Ann Microbiol. 2016;66(1):151–64.
Yang W, Wang L, Hu Q, Pei F, Mugambi MA. Identification of bacterial composition in freeze-dried Agaricus bisporus during storage and the resultant odor deterioration. Front Microbiol. 2019;10(FEB):1–12.
Pecchia J, Cortese R, Albert I. Investigation into the microbial community changes that occur in the casing layer during cropping of the white button mushroom, Agaricus bisporus. In: Singh M, editor. Proceedings of the 8th International Conference on Mushroom Biology and Mushroom Products, New Delhi; 2014. p. 309–13.
McGee CF, Byrne H, Irvine A, Wilson J. Diversity and dynamics of the DNA-and cDNA-derived compost fungal communities throughout the commercial cultivation process for Agaricus bisporus. Mycologia. 2017;109(3):475–84.
Braat N, Koster MC, Wösten HAB. Beneficial interactions between bacteria and edible mushrooms. Fungal Biol Rev. 2022;39:60–72.
Sigler WV, Turco RF. Erratum: The impact of chlorothalonil application on soil bacterial and fungal populations as assessed by denaturing gradient gel electrophoresis (Applied Soil Ecology (2002) 21 (107-118) PII: S0929139302000884). Appl Soil Ecol. 2003;22(2):193.
Nishiyama T, Ueki A, Kaku N, Watanabe K, Ueki K. Bacteroides graminisolvens sp. nov., a xylanolytic anaerobe isolated from a methanogenic reactor treating cattle waste. Int J Syst Evol Microbiol. 2009;59(8):1901–7.
Carrasco J, García-Delgado C, Lavega R, Tello ML, De Toro M, Barba-Vicente V, et al. Holistic assessment of the microbiome dynamics in the substrates used for commercial champignon (Agaricus bisporus) cultivation. Microb Biotechnol. 2020;13(6):1933–47.
Gandy DO, Spencer DM. Fungicides for the control of Mycogone perniciosa (Magn.), the cause of wet bubble on the cultivated mushroom. Sci Hortic. 1978;8:307–13.
Taparia T, Hendrix E, Nijhuis E, de Boer W, van der Wolf J. Circular alternatives to peat in growing media: a microbiome perspective. J Clean Prod. 2021;327(October):129375.
Zhang X, Zhong Y, Yang S, Zhang W, Xu M, Ma A, et al. Diversity and dynamics of the microbial community on decomposing wheat straw during mushroom compost production. Bioresour Technol. 2014;170:183–95.
Rocchi S, Valot B, Reboux G, Millon L. DNA metabarcoding to assess indoor fungal communities: electrostatic dust collectors and Illumina sequencing. J Microbiol Methods. 2017;139(January):107–12.
Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013;41(1):1–11.
White TJ, Bruns T, Lee S, Taylor J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. PCR Protoc. 1990;1:315–22.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13(7):581.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high- throughput community sequencing data intensity normalization improves color calling in SOLiD sequencing. Nat Publ Gr. 2010;7(5):335–6.
Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37(8):852–7.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41(Database issue):D590–6.
Abarenkov K, Zirk A, Piirmann T, Pöhönen R, Ivanov F, Nilsson RH, et al. UNITE QIIME release for Fungi [data set]: UNITE Community; 2020. https://doi.org/10.15156/BIO/786385.
Mandal S, Van Treuren W, White RA, Eggesbø M, Knight R, Peddada SD. Analysis of composition of microbiomes: a novel method for studying microbial composition. Microb Ecol Heal Dis. 2015;26(0):1–7.
Faith DP. Conservation evaluation and phylogenetic diversity. Biol Conserv. 1992;61(1):1–10.
Pielou EC. The measurement of diversity in different types of biological collections. J Theor Biol. 1966;13(C):131–44.
Shannon CE, Weaver W. The mathematical theory of communications. Urbana: University of Illinois Press; 1949.
Kruskal WH, Wallis WA. Use of ranks in one-criterion variance analysis. J Am Stat Assoc. 1952;47(260):583–621.
Weiss S, Xu ZZ, Peddada S, Amir A, Bittinger K, Gonzalez A, et al. Normalization and microbial differential abundance strategies depend upon data characteristics. Microbiome. 2017;5(1):27.
Sorenson T. A method of establishing groups of equal amplitude in plant sociology based on similarity of species content. Kongelige Danske Videnskabernes Selskab. 1948;5(1–34):4–7.
Lozupone C, Knight R. UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol. 2005;71(12):8228–35.
RStudio Team. RStudio: integrated development for R. Boston: RStudio, PBC; 2020. http://www.rstudio.com/
Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, et al. Welcome to the Tidyverse. J Open Source Softw. 2019;4(43):1686.
Bisanz JE. Qiime2R: importing QIIME2 artifacts and associated data into R sessions. v0.99. 2018. https://github.com/jbisanz/qiime2R.
This work was co-funded by the FEDER Operational Program La Rioja 2014–2020 (Project CT21_04), the Department of Economic Development and Innovation of La Rioja. Res. n° 1168/2018 and the H2020 Marie Sklodowska-Curie Actions, Grant agreement no. 742966.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
ANCOM Statistical analysis on the bacteria population with the significant genus data-set DAYS.
ANCOM Statistical analysis on the fungi population with the significant genus data-set DAYS.
ANCOM Statistical analysis on the bacteria population with the significant genus data-set TREATMENTS-D16.
Diversity significant boxplot for 16S analysis by DAY. Right to left: D00, D04, D07 and D16. With Faith Phylodiversity (1A), Pielou’s species evenness (1B) and Shannon’s index (1C). Three-dimensional PCoA visualized using Emperor, built using the Bray Curtis distance matrix (1D) and the unweighted UniFrac distance matrix (1E).
Diversity significant boxplot for 16S analysis by TREATMENT. Right to left: Chlorothalonil, Control and Metrafenone. With Faith Phylodiversity (2A), Pielou’s species evenness (2B) and Shannon’s index (2C). Three-dimensional PCoA visualized using Emperor, built using the Bray Curtis distance matrix (2D) and the unweighted UniFrac distance matrix (2E).
Diversity significant boxplot for 16S analysis by TREATMENT. Right to left: Chlorothalonil, Control and Metrafenone. With Faith Phylodiversity (3A), Pielou’s species evenness (3B) and Shannon’s index (3C). Three-dimensional PCoA visualized using Emperor, built using the Bray Curtis distance matrix (3D) and the unweighted UniFrac distance matrix (3E).
Diversity significant boxplot for ITS2 analysis by DAY. Right to left: D00, D04, D07 and D16. With Faith Phylodiversity (4A), Pielou’s species evenness (4B) and Shannon’s index (4C). Three-dimensional PCoA visualized using Emperor, built using the Bray Curtis distance matrix (4D) and the unweighted UniFrac distance matrix (4E).
Diversity significant boxplot for ITS2 analysis by TREATMENT. Right to left: Chlorothalonil, Control and Metrafenone. With Faith Phylodiversity (5A), Pielou’s species evenness (5B) and Shannon’s index (5C). Three-dimensional PCoA visualized using Emperor, built using the Bray Curtis distance matrix (5D) and the unweighted UniFrac distance matrix (5E).
Diversity significant boxplot for ITS2 analysis by TREATMENT. Right to left: Chlorothalonil, Control and Metrafenone. With Faith Phylodiversity (6A), Pielou’s species evenness (6B) and Shannon’s index (6C). Three-dimensional PCoA visualized using Emperor, built using the Bray Curtis distance matrix (6D) and the unweighted UniFrac distance matrix (6E).
Rarefaction plots of all 23 samples for 16S analysis (sample 26 excluded). 7A: Rarefaction curves (Shannon’s index on Y axis), 7B: Rarefaction curves (Number of observed features on Y axis).
Rarefaction plots of all 24 samples for ITS2 analysis. 8A: Rarefaction curves (Shannon’s index on Y axis), 8B: Rarefaction curves (Number of observed features on Y axis).
Heatmap with the bacterial taxonomy at Species level of data-set DAYS.
Heatmap with the bacterial taxonomy at Species level of data-set TREATMENTS D07.
Heatmap with the bacterial taxonomy at Species level of data-set TREATMENTS D16.
Heatmap with the fungal taxonomy at Species level of data-set DAYS.
Heatmap with the fungal taxonomy at Species level of data-set TREATMENTS D07.
Heatmap with the fungal taxonomy at Species level of data-set TREATMENTS D16.
About this article
Cite this article
Tello Martín, M.L., Lavega, R., Carrasco, J.C. et al. Influence of Agaricus bisporus establishment and fungicidal treatments on casing soil metataxonomy during mushroom cultivation. BMC Genomics 23, 442 (2022). https://doi.org/10.1186/s12864-022-08638-x