Transcriptome Analysis and Expression Profiling of Molecular Responses to Cd Toxicity in Morchella spongiola

Abstract

Morchella is a genus of fungi with the ability to concentrate Cd both in the fruit-body and mycelium. However, the molecular mechanisms conferring resistance to Cd stress in Morchella are unknown. Here, RNA-based transcriptomic sequencing was used to identify the genes and pathways involved in Cd tolerance in Morchella spongiola. 7444 differentially expressed genes (DEGs) were identified by cultivating M. spongiola in media containing 0.15, 0.90, or 1.50 mg/L Cd²⁺. The DEGs were divided into six sub-clusters based on their global expression profiles. GO enrichment analysis indicated that numerous DEGs were associated with catalytic activity, cell cycle control, and the ribosome. KEGG enrichment analysis showed that the main pathways under Cd stress were MAPK signaling, oxidative phosphorylation, pyruvate metabolism, and propanoate metabolism. In addition, several DEGs encoding ion transporters, enzymatic/non-enzymatic antioxidants, and transcription factors were identified. Based on these results, a preliminary gene regulatory network was firstly proposed to illustrate the molecular mechanisms of Cd detoxification in M. spongiola. These results provide valuable insights into the Cd tolerance mechanism of M. spongiola and constitute a robust foundation for further studies on detoxification mechanisms in macrofungi that could potentially lead to the development of new and improved fungal bioremediation strategies.

Keywords:

• Cadmium stress
• transporter protein
• enzymatic antioxidant
• non-enzymatic antioxidant
• gene regulatory network

1. Introduction

Cadmium (Cd) is one of five heavy metals with long biological half-lives that are commonly found as environmental pollutants [1]. Exposure to Cd causes a variety of problems in microorganisms including reduced accumulation of microbial biomass carbon (C_mic), changes in microbial community structure and losses of diversity, and sharp reductions in enzyme activity [2]. In recent years, several studies have shown that both cultured and wild fruiting bodies of mushrooms such as Lentinus edodes [3], Agaricus subrufescens [4], Galerina vittiformis [5], and Agaricus bisporus [6] can absorb Cd from polluted soil. To date this effect has only been studied in the environment with a single Cd content from 0.366 mg/L to 52.9 mg/kg [7,8]. There have been few detailed studies on the mechanisms of Cd uptake and resistance in microorganisms under different doses of Cd exposing. Understanding how edible fungi resist Cd stresses could clarify fungal heavy metal detoxification mechanisms and enable the development of more efficient bioremediation strategies for treating soil with severe heavy metal contamination.

Known molecular mechanisms of Cd stress generally involve changes in the expression of genes involved with DNA damage repair, inhibition of apoptosis and autophagy, oxidative stress management, and interaction with biological elements [9]. For instance, in Volvariella volvacea, Cd stress induces selective overexpression of genes encoding metallothioneins (MTs) and heat shock proteins (HSPs), along with genes involved in oxidative stress responses [7]. Cd resistance in Bacillus megaterium NCT-2 and Bacillus paranthracis NT1 may be regulated by the genes, like cada, zitB, khtT, bshA, trkA, czcD and abshA [10]. Additionally, some proteins related to DNA synthesis and repair, including ribosomal proteins, phosphoglycerate kinase, and inorganic pyrophosphatase were reportedly up-regulated by Cd stress [11]. Cd concentration of the fungal growth environment is usually very low and can reach the level of biotoxicity only in the vicinity of natural or man-made point sources. How these subtoxic concentrations of Cd affect the transcriptional level of genes remains to be studied in more detail.

Morchella species are known worldwide for their flavor and diverse bioactivities and were found to be enriched in Cd when cultivated in Cd-contaminated soil [12–15]. Cd reportedly stimulated mycelial growth in Morchella esculenta at concentrations of 0–2 mg/L [16]. However, concentrations above 50 mg/L strongly inhibited mycelial growth in Morchella importuna [16,17]. Preliminary experimental studies showed that the presence of Cd (>0.9 mg/L) suppressed the growth of M. spongiola in culture medium and induced changes in hyphal micro-morphology, including twisting, folding, and kink changing. To clarify the molecular mechanisms by which Cd suppresses growth in M. spongiola, mycelia were cultivated in media with different concentrations of CdCl₂ and then analyzed using transcriptome RNA-seq sequencing to identify gene expression patterns associated with different Cd concentrations. The results obtained provide new insights into the mechanisms of Cd tolerance in M. spongiola.

2. Materials and methods

2.1. Mycelium culture and Cd treatment

Morchella spongiola M12-10 was obtained from the Laboratory of Microbiology, College of Eco-Environmental Engineering, Qinghai University, China. Mycelium samples were cultured in Potato Dextrose Agar medium (glucose 20.0 g/L, potato extraction 200.0 g/L, and agar 17.0 g/L) in an incubator at 20 °C and 50% humidity for 5 days. Three Cd concentrations (0.15, 0.90, and 1.50 mg/L CdCl₂) and a control treatment without Cd were used to stress the mycelia of M. spongiola. Media containing these different concentrations of Cd were prepared as described previously [17]. Ten replicate plates were prepared and analyzed for each Cd treatment including the control group. After 5 days’ Cd exposure, the mycelia were collected, mixed and pooled in Eppendorf-tubes (1.5 mL), rapidly frozen using liquid nitrogen, and stored at −80 °C for RNA extraction.

2.2. RNA extraction and illumina sequencing

For each of Cd treatments and the control, RNA was extracted from 20 mg of mixed mycelia. Total RNA was extracted from each sample using Total RNA Extractor (Trizol) (Sangon Biotech, China) following the manufacturer's protocol. After assessed the concentration and quality of total RNA by Invitrogen Qubit 2.0 Fluorometer (Thermo Fisher Scientific, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA), the VAHTS™ mRNA-seq V2 Library Prep Kit for Illumina kit (Vazyme Biotech, Nanjing, China) was used to RNA fragmentation, double-stranded cDNA synthesis, repair sticky ends, terminal dA-Tailing, joint connection, purification of ligated products, sorting of fragment size and the library amplification. Then, cDNA was subjected to PCR amplification based on the size of the selected fragment. Finally, we performed Illumina double-ended sequencing (150 bp per read) on an Illumina Hiseq X^Ten instrument after accurately quantifying the recovered DNA (Qubit 2.0 DNA detection kit; Life Biotech, Shanghai, China). The original image data file obtained by Illumina Hiseq recognition was transformed into the original sequencing sequence (Sequenced Reads) base on the CASAVA base recognition (Base Calling) analysis. Then the phred quality score (Q) was used to map the read. The mapping relationship is Q = −10^{log (P)}, while P was the base-calling error probabilities. The base-calling error probabilities were 90%, 99%, 99.9%, and 99.99% when the phred quality scores were 10, 20, 30, and 40, respectively.

2.3. Read filtering and denovo assembly

The raw sequence data was dominated by linkers and low-quality sequences. Therefore, the quality of the sequenced raw data was evaluated using FastQC, trimmed using Trimmomatic, mixed and spliced using Trinity to obtain clean data. Trinity was then used to denovo assemble the clean data into a transcript. The resulting spliced transcript was used as a reference sequence for redundant assessment using RSeQC, and the longest transcript in each transcript cluster was treated as a unigene in subsequent analyses.

2.4. Differential gene expression analysis and functional annotation

Gene expression counts were generated using Salmon. To compare the expression of different genes, we expressed the relative abundance of each transcript in an RNA pool in units of Transcripts Per Million (TPM). Read counts were normalized using TMM and differentially expressed genes (DEGs) were identified with DEGseq, using a cutoff threshold of q-value ≦ 0.05 and |FoldChange| ≧ 2 to determine significant differences in expression. The annotation information included a cluster of orthologous group (COG) functional annotations, protein functional annotations, and gene ontology (GO) functional annotations for unigenes based on BLASTX alignments to the NR, Swiss-Prot, KEGG and COG protein databases. GO terms were obtained using the ClusterProfiler software. Additionally, pathway analysis was performed using the KEGG mapping method. The unigene sequences were mapped to KEGG biochemical pathways based on their EC numbers using the KEGG Automatic Annotation Server.

2.5. Quantitative real-time PCR validation

To verify the reliability of the DEGs identified from the RNA-Seq data, we performed qRT-PCR on the same samples used for transcriptome profiling. Twenty DEGs were randomly selected and primers were designed for them using Primer 5.0. The 18S rRNA gene was used as the reference gene, and the primers used for qRT-PCR are listed in a supplementary table (Table S1). Then, qRT-PCR was performed using a Bio-rad CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, USA) and 2*Taq qPCR Master Mix (Sangon Biotech, Shanghai, China). Reactions were performed with three technological and three biological replicates using 1.0 µL of each primer, 9.5 µL ddH₂O, 1.0 µL of cDNA, and 12.5 µL of 2*Taq qPCR Master Mix at a final volume of 25 µL. The PCR program was as follows: 1 min at 94 °C followed by 35 cycles of amplification at 95 °C for 15 s, then 59 °C for 15 s, and 72 °C for 1 min, with a final 5 min period at 72 °C. Fold changes in gene expression were calculated using the 2^−△△CT method, and a correlation analysis was performed to study the relationships between the expression strength of the DEGs as determined by qRT-PCR and RNA-Seq [18].

2.6. Determination of SOD and H₂O₂

The collected mycelia were homogenized in PBS (50 mm pH 7.8) at the ratio of 1:9 (biomass: volume), and then centrifuged at 10,000 rpm for 10 min at 4 °C. The supernatant was obtained as the extract of M. spongiola mycelia. SOD activity and H₂O₂ accumulation in the extract were determined according to the assay kits of SOD and H₂O₂ (Nanjing Jiancheng Bioengineering Institute, Nanjing, China).

3. Results and discussion

3.1. Quality and reliability of transcriptome sequencing and unigene assembly

To investigate gene expression changes in M. spongiola under Cd stress, RNA libraries corresponding to the control and 0.15 mg/L, 0.90 mg/L, and 1.50 mg/L Cd treatments were sequenced using the Illumina Hiseq platform, and 30,653,804, 71,591,348, 127,771,204, 77,562,142 clean reads for each treatment were obtained (Table S2). Over 98.64% of these reads were clean at the Q20 level and over 95.05% at the Q30 level, indicating that the sequencing data was of sufficient quality to justify further analysis. Additionally, the GC content of the sequence data was 51%, indicating that the sequencing data were stable (Table S2). The results of a pollution assessment of the sequencing data in Table S3 indicated that sequences from other species were not represented in the data. Overall, these results demonstrate the credibility of the transcriptome sequencing data.

3.2. Differential gene expression

A total of 7444 significant DEGs were identified, of which 564, 1023, and 1862 were significantly up-regulated and 951, 2514, and 530 were significantly down-regulated under the 0.15 mg/L, 0.90 mg/L, and 1.50 mg/L Cd treatments, respectively (q ≦ 0.05) (Figure 1(A)). 20.35%, 47.51%, and 32.13% DEGs under above three Cd concentrations, demonstrated that M. spongiola strongly alters its gene expression to respond to the toxicity stress by different concentration Cd exposing. In addition, 305 genes were co-differentially expressed under all three Cd treatments, and they are likely to play key roles in resisting damage caused by Cd stress (Table S4). The set of Cd-related coregulatory genes was enriched in the functions “integral membrane component,” “nucleusin,” “oxidoreductase activity,” “ATP binding,” “transferase activity,” “carbohydrate metabolic process,” “response to stress” and “transmembrane transport,” suggesting that M. spongiola resists Cd toxicity by regulating the expression of genes that affect cell membrane composition and transport, stress responses, and redox metabolism.

Figure 1. Identification and classification of DEGs under the 0.15, 0.90, 1.50 mg/L Cd treatments in Morchella spongiola. (A) The number of up- and down-regulated DEGs for each treatment. (B) Venn diagram showing the total number of DEGs annotated with the function “response to cadmium ion” (GO: 0046686) for the three treatments.

3.3. Functional enrichment analysis of DEGs

Table 1 shows the 20 most strongly enriched GO annotations for the DEGs. For the 0.15 mg/L treatment, the most strongly enriched gene function categories were “catalytic activity,” “single-organism metabolic process,” “intrinsic component of membrane,” and “integral component of membrane.” However, cellular component annotations were more strongly enriched for the 0.90 mg/L treatment, while the most strongly enriched terms for the 1.50 mg/L Cd treatment were related to biological processes. These results suggest that M. spongiola may prevent exogenous low concentration Cd from entering the cell by promoting enzyme-catalyzed reactions and metabolic processes that maintain membrane stability and integrity. It should be noted that there is a specific GO term (GO: 0046686) for genes associated with Cd stress responses, and there were 4 such DEGs for the 0.15 mg/L treatment, 23 for the 0.90 mg/L treatment, and 20 for the 1.50 mg/L treatment (Figure 1(B)). One of these DEGs was up-regulated under all three Cd stress treatments. This gene, TRINITY_DN54495_c5_g6, encodes a thaumatin-like protein; genes of this type are known to be involved in host defense in fungi, plants, and animals [19].

Table 1. The 20 most significantly enriched GO terms associated with DEGs for each of the three Cd stress treatments.

Treatment group	Term description	Term_type	Up-regulated DEGs number	Down-regulated DEGs number
0.15 mg/L	Catalytic activity	Molecular function	338	193
	Single-organism metabolic process	Biological process	182	92
	Intrinsic component of membrane	Cellular component	137	113
	Integral component of membrane	Cellular component	129	111
	Cytosol	Cellular component	147	66
	Hydrolase activity	Molecular function	128	65
	Small molecule metabolic process	Biological process	141	41
	Organonitrogen compound metabolic process	Biological process	133	40
	Oxidoreductase activity	Molecular function	80	56
	Single-organism biosynthetic process	Biological process	92	38
	Organonitrogen compound biosynthetic process	Biological process	96	24
	Organic substance catabolic process	Biological process	83	36
	Organic acid metabolic process	Biological process	88	24
	Oxoacid metabolic process	Biological process	88	23
	Carboxylic acid metabolic process	Biological process	86	21
	Carbohydrate metabolic process	Biological process	65	28
	Organophosphate metabolic process	Biological process	46	23
	Transmembrane transport	Biological process	40	27
	Single-organism catabolic process	Biological process	45	21
	Small molecule biosynthetic process	Biological process	52	10
0.90 mg/L	Cell	Cellular component	1403	727
	Cell part	Cellular component	1401	725
	Cellular process	Biological process	1219	659
	Intracellular	Cellular component	1337	700
	Intracellular part	Cellular component	1319	697
	Organelle	Cellular component	1129	665
	Cytoplasm	Cellular component	1097	672
	Intracellular organelle	Cellular component	1119	633
	Metabolic process	Biological process	1022	544
	Organic substance metabolic process	Biological process	968	511
	Binding	Molecular function	963	521
	Cellular metabolic process	Biological process	940	519
	Membrane-bounded organelle	Cellular component	1074	518
	Primary metabolic process	Biological process	892	493
	Intracellular membrane-bounded organelle	Cellular component	1044	459
	Single-organism process	Biological process	908	428
	Cytoplasmic part	Cellular component	908	572
	Catalytic activity	Molecular function	882	246
	Single-organism cellular process	Biological process	722	359
	Nitrogen compound metabolic process	Biological process	662	393
1.50 mg/L	Metabolic process	Biological process	169	908
	Cytoplasm	Cellular component	155	976
	Organic substance metabolic process	Biological process	151	850
	Cellular metabolic process	Biological process	146	823
	Primary metabolic process	Biological process	129	801
	Cytoplasmic part	Cellular component	130	837
	Nitrogen compound metabolic process	Biological process	96	591
	Biosynthetic process	Biological process	81	552
	Organic substance biosynthetic process	Biological process	74	545
	Cellular biosynthetic process	Biological process	76	519
	Macromolecular complex	Cellular component	45	490
	Organonitrogen compound metabolic process	Biological process	53	413
	Cellular nitrogen compound biosynthetic process	Biological process	42	414
	Protein metabolic process	Biological process	23	431
	Single-organism metabolic process	Biological process	93	354
	Cellular protein metabolic process	Biological process	22	412
	Cytosol	Cellular component	55	377
	Organonitrogen compound biosynthetic process	Biological process	32	323
	Response to stress	Biological process	20	303
	Cellular amide metabolic process	Biological process	11	265

Pathway reconstruction analysis revealed that DEGs identified for each Cd stress treatment were enriched in several metabolic pathways, including “Glycolysis/Gluconeogenesis,” “Pentose phosphate pathway,” “Pentose and glucuronate,” “Glycerolipid metabolism,” “Styrene degradation,” and “Biosynthesis of amino acids” (Figure 2). Pathways enriched in the 0.15 mg/L treatment included carbon metabolism, starch and sucrose metabolism, and aminoacyl-tRNA biosynthesis. Pathways under the 0.90 mg/L treatment enriched in MAPK signaling, carbon metabolism, and the cell cycle, while it include valine, leucine and isoleucine degradation, fatty acid metabolism, and lysine biosynthesis under the 1.50 mg/L treatment.

Figure 2. Significant pathways and numbers of DEGs in different groups. Significantly enriched pathways and numbers of DEGs associated with those pathways under the (A) 0.15 mg/L, (B) 0.90 mg/L treatment, and (C) 1.50 mg/L Cd treatments.

3.4. Expression profile of DEGs

To shed further light on the changes in gene expression induced by Cd stress in M. spongiola, the Gplot package was used to cluster the expression patterns of significant DEGs. The 7444 DEGs were divided into 6 different subclusters based on their expression patterns (Figure 3). The mRNA-level expression of genes in subcluster 1 fell gradually as the Cd concentration increased; genes in this subcluster were expressed most strongly under the control treatment and least strongly at 1.50 mg/L Cd. DEGs in subcluster 2 were generally up-regulated under the control, 0.15 mg/L, and 0.90 mg/L treatments, but significantly down-regulated at 1.50 mg/L Cd. DEGs in subcluster 3 were down-regulated under the control and 0.15 mg/L treatments but up-regulated under the 0.90 mg/L and 1.50 mg/L treatments. DEGs in subcluster 4 were downregulated under the control, 0.15 mg/L, and 0.90 mg/L treatments, but strongly up-regulated under the 1.50 mg/L treatment. DEGs in subcluster 5 were up-regulated under the control and 0.15 mg/L treatments but down-regulated under the 0.90 mg/L and 1.50 mg/L treatments. The mRNA expression of DEGs in subcluster 6 generally increased with the Cd concentrations. GO enrichment and KEGG pathway analysis showed that genes in subcluster 1 involved in biosynthesis of amino acids, starch and sucrose metabolism, aminoacyl-tRNA biosynthesis, and proteasome (Figure S1), in subcluster 2 involved in RNA processes and MAPK signaling pathway (Figure S2), in subcluster 3 and 4 involved in protein biosynthesis in the ribosome (Figures S3 and S4), in subcluster 5 involved in amino metabolism (Figure S5), in subscluster 6 involved in oxidative phosphorylation (Figure S6).

Figure 3. Heat-map showing changes in the expression of DEGs in six sub-clusters at increasing Cd concentrations.

3.5. Validation of transcriptome analysis

To verify the reliability of the RNA-seq data, intracellular SOD activity and H₂O₂ accumulation were determined, 20 RNA samples were also randomly selected and their expression profiles in M. spongiola under Cd stress were investigated by qRT-PCR. It was clear that SOD activity of the Cd-treated samples were higher than those of the control (Figure 4(A)). The content of H₂O₂ in 0.15 mg/L, 0.90 mg/L and 1.50 mg/L Cd-treated group was 3.91, 5.40 and 8.75 times higher than that in control group, respectively (Figure 4(B)). Furthermore, the RNA-seq and qRT-PCR expression data were strongly correlated (Figure 4(C); Table S5); the Pearson’s correlation coefficient for the two datasets (r = 0.9797) indicated good concordance between the two techniques. We can therefore conclude that the RNA-seq data accurately reflect the transcriptomic changes in M. spongiola induced by Cd stress.

Figure 4. Validation of SOD, H₂O₂ and qRT-PCR. (A) Intracellular SOD activity of Morchella spongiola exposed to 0.15, 0.90, 1.5 mg/L Cd; The inserted figure showed the significant correlation between initial Cd concentration and SOD in M. spongiola. (B) Intracellular H₂O₂ accumulation of M. spongiola exposed to 0.15, 0.90, 1.5 mg/L Cd; The inserted figure showed the significant correlation between initial Cd concentration and H₂O₂ in M. spongiola. (C) Correlation between the RNA-seq and qRT-PCR data. Each point represents the log₂FC values of one of the 20 chosen DEGs based on its relative expression as determined by RNA-seq and qRT-PCR.

3.6. Expression of transporter proteins in M. spongiola under Cd stress

Heavy metal transport proteins play a key role in microbiological metabolism under Cd stress. In particular, ATP binding cassette transporters (ABC transporters) responsible for maintaining metal ion homeostasis by transporting metal ions across the cell membrane, either as naked ions or in the form of various complexes [20,21]. In this work, 7 ABC transporters were identified in M. spongiola; their mean change in expression was 1.006< log₂FC <15.776 (Table S6). Notably, the ABC1 family proteins YPL109C and MRP CDR2 were up-regulated under the 0.15 mg/L Cd treatment, suggesting that these two genes may play key roles in resistance to low exogenous Cd concentrations in M. spongiola. ABC transporter F family member 5 was also strongly up-regulated under the 1.50 mg/L treatment (log₂FC = 15.776), indicating that ABCG5 could play an important role in enhancing resistance to Cd toxicity in M. spongiola. Copper transport protein (CTR2) is mainly found in intracellular membranes, where it acts as a vacuolar/vesicular transporter responsible for regulating the intracellular copper concentration by moving copper from vesicular organelles to the cytoplasm [22]. Three CTR-encoding genes, were also upregulated under the 0.90 mg/L treatment, and are probably responsible for regulating intracellular copper homeostasis under Cd stress. However, there was no appreciable transcript-level expression of genes encoding CTR proteins under the 1.50 mg/L Cd treatment, suggesting that severe Cd stress may abolish intracellular copper homeostasis or activate other response mechanisms (Table S6). Additionally, 3 Zinc-regulated transporters (ZRT) genes were up-regulated under the 0.15 mg/L treatment, suggesting that low concentrations of Cd may cause competitive zinc deficiency, provoking the expression of zinc-regulated transporter genes to maintain intracellular zinc homeostasis (Table S6).

3.7. Expression of enzymatic and non-enzymatic antioxidant genes in M. spongiola under Cd stress

Enzymes play important roles in the abiotic stress responses of microorganisms. Series of up-regulated enzymes genes, such as sodA, pod, cat, Fsd 1, Msd 3, Apx 1, and Epgst, were reportedly increased after exposure to Cd²⁺ [23–25]. In our study, three CuSod genes were strongly up-regulated, with log₂FC values of 15.99–17.46 (Table S7). Conversely, the upregulation of the MnSod gene was less pronounced. This suggests that CuSod may play a particularly important role in M. spongiola under Cd stress. Two cat genes were up-regulated under the 0.90 mg/L treatment, indicating that cat genes may only contribute to ROS scavenging at this higher level of Cd stress. The pod gene was only up-regulated under the 1.50 mg/L treatment, suggesting that it may only contribute to ROS scavenging under especially severe Cd stress. Additionally, four gst genes and two gsh-px genes were significantly up-regulated under Cd stress, but gst1 was down-regulated. It thus appears that gst genes contribute significantly to Cd stress resistance in M. spongiola. Based on the expression difference of antioxidant genes and physilogical changes of SOD, H₂O₂ (Figure 4(A,B)), we hypothesize that the response to low levels of Cd stress in M. spongiola depends primarily on the heavy metal transport system, while ROS scavenging by the antioxidant system plays an increasingly important role as the severity of the Cd stress increases.

Non-enzymatic antioxidant systems are also important in abiotic stress responses. In fungi, HSPs promote the folding, stabilization, transport and degradation of proteins. As such, they are involved in regulating cell cycle progression and the activation of many key signal transduction enzymes [26]. All three Cd-treatments applied in this work significantly altered the expression of multiple genes, including 17 HSP, 8 TRX, and 5 MT genes (Table S8). Under 0.15 mg/L Cd treatment, five hsp genes and two trx genes, including hsp12, hsp70, hsp98, hspssc1, and trx were upregulated in M. spongiola. Under the 0.90 mg/L treatment, 8 HSP genes including hsp12, hsp90, and LMW-hsp were differentially expressed, with the highest log₂FC being 19.96, indicating a strong response to Cd toxicity. Thus, we speculate that hsp70 and hsp12 may be involved in cellular morphogenesis and cell adhesion in M. spongiola under Cd stress, while hsp90 may be important for maintaining genomic stability, and several LMW HSPs may be primarily responsible for preventing irreversible protein polymerization. Five trx genes were also strongly up-regulated, with log₂FC values of 5.03–16.67, and five mt genes were very strongly up-regulated, with log₂FC values of 15.99- 19.60 under both the 0.90 mg/L and 1.50 mg/L treatments. This suggests that trx and mt genes may increase the Cd tolerance of M. spongiola under severe Cd stress by altering the redox state of cells and chelating ions.

3.8. Expression of transcription factors in M. spongiola under Cd stress

Cd-responsive transcription factors have been shown to participate in Cd stress responses. Several transcription factor genes exhibited significant changes in M. spongiola under Cd stress, including gata, atf21, btf3, myb, and sfc1 (Table S9). In this work, two TF genes (namely atf21 and myb), were up-regulated under the 0.15 mg/L treatment, while btf3 was up-regulated under the 0.90 mg/L treatment and a myb44-like gene was up-regulated under the 1.50 mg/L Cd treatment. It was therefore speculated that these transcription factors may responsible for Cd stress adaptation in M. spongiola.

3.9. Mechanisms of Cd stress resistance in M. spongiola

To clarify the mechanisms of resistance to heavy metal toxicity in macro-fungi, we conducted a transcriptomic study on concentration-dependent responses to three levels of Cd stress in M. spongiola. First of all, heavy metal transport proteins play a key role in microbiological metabolism under Cd stress. In this study, a large number of transporter genes were remarkably induced in different Cd-treatment groups, including 7 ATP binding cassette transporter (ABC transporter), 3 Multidrug resistance protein (MRP), 4 Copper transporter (CTR), 4 Zinc transporter (ZRT), 3 Heavy metal transporter (HMT), 1 Cation transport ATPase (CE) and 1 Iron transport multicopper (ILT) (Table S6). ABC transporters, as one of the largest families of microbial membrane proteins, shows great important in resisting virulence and drug, especially in the function of maintaining metal ion homeostasis by transporting metal ions across the cell membrane in ionic/complex forms [20,21]. It was reported that the expression of 3 ABC genes (EpABC2.1, EpABC3.1, or EpABC4.1) of dominant dark septate endophytes (DSEs) could promote the restoring growth of metal-sensitive mutant by alleviating the inhibitory activity in yeast caused by Pb, Zn and Cd, and further altered toxic ions uptake and accumulation [27]. A specific ABC-transporter of Escherichia coli could mediate the entry of the Fe²⁺ bounding siderophores through the inner membrane of the cell from the outer membrane under the help of FepB and FhuD chaperones [28]. Multidrug resistance protein (MRP/ABCC) is a broad-spectrum ABC transporter that plays a major role in defense against environmental toxicants, in addition to contributing toward multidrug resistance of certain types of malignancy [29]. The ABC transporters CDR2 was reported to play a key role in azole resistance of yeast Candida [30]. In our study, ABC1 family protein YPL109C (DN52619_c4_g1) and MRP CDR2 (DN46433_c0_g1) were up-regulated at 0.15 mg/L Cd-treatment, suggesting that these genes might play a key role in low concentration Cd resistance in M. spongiola. One ABC gene (DN53140_c5_g2) was significantly increased at a growing rate from 0.90 mg/L to 1.50 mg/L, indicating their strong response to Cd stress under the higher concentrations. ABC transporter G family member 5 (DN50046_c0_g3) has been found to be highly upregulated under 1.50 mg/L Cd stress, indicating that ABCG5 might play an important role in enhancing the resistance to Cd toxicity in M. spongiola. Copper transport protein (CTR2) is mainly found in intracellular membranes as a vacuolar/vesicular transporter, which responsible for regulating the intracellular copper concentration by moving copper from vesicular organelles to the cytoplasm [22]. CTR2 of Saccharomyces cerevisiae has been proved to regulate intercellular Cu²⁺ homeostasis under low copper conditions [31]. CTR4, as one of the high-affinity Cu²⁺ transporter, mainly responsible for the intercullar Cu²⁺ homeostasis and virulence in Cryptococcus neoformans [32]. ATOX1, as one of metallochaperone in organisms, directly interacts with the Cu²⁺-transporting ATPases and plays a foundational role in Cu²⁺ homeostasis [33,34]. Kim et al. reported that ATOX1 is a necessary chaperon for MAPK signaling and NF-κB pathways in BRAF mutation-positive melanoma [35,36]. In our study, CTR ATOX1 (DN52132_c3_g5) and CTR (DN3415_c0_g1) were up-regulated under 0.90 m/L Cd stress, suggesting that they probably responsible for regulating intracellular Cu homeostasis to defense the cellular damage caused by Cd stress. While the gene transcripts of encoding CTR did not show significant expression in 1.50 mg/L Cd-treatment group, indicating that intracellular copper homeostasis might have been destroyed by ultra-high concentration Cd stress or other response mechanisms have been induced by ultra-high concentration Cd stress. Zinc-regulated transporters (ZRT), as a member of the ZIP family, are responsible for the uptake and transport of divalent metal ions to maintain ion homeostasis in the organisms [37]. In microorganisms, zinc transporter expression is mainly regulated to maintain intracellular zinc homeostasis [38]. 2 up-regulated ZRT gene (F751_0752 and F751_4157) were predicted resist Cd stress at the early growth stage of microalgae [19]. The levels of the ZRT1, ZRT3, and FET4 metal transporters of Saccharomyces cerevisiae are increased during zinc deficiency by Zap1 to facilitate zinc uptake and mobilize zinc stored in the vacuole [39]. In our study, 3 ZRT genes were up-regulated under 0.15 mg/L treatment, suggesting that low concentrations of Cd might lead to competitive zinc deficiency, then induces expression of zinc-regulated transporter genes, which helped to maintain intracellular zinc homeostasis. Lu et al. reported that both calcium ATPase (Ca-ATPase) and cation efflux ATPase (CE) plays a great role in moving Cd into cells [18]. In our study, DN54570_c1_gi and DN45045_c0_g1 genes were upregulated both in 1.5 mg/L treatment and 1.5 mg/L group, indicating they might have played an important role in absorption and transport of Cd in M. spongiola. Taken together, the transcripts difference in the expression of transporter genes revealed that M. spongiola has different coping mechanisms in different concentrations of Cd-stress.

Secondly, enzymatic antioxidant genes in M. spongiola play an important role in abiotic stress. Cd is a common cause of reactive oxygen species (ROS) outburst in organisms, which induces severe damage to cell structures and ultimately results in cell death [40]. To scavenge ROS, organisms have been developed complex antioxidant systems that include both non-enzymatic and enzymatic antioxidants [41]. The antioxidant enzymes such as superoxide dismutase (SOD), peroxidase (POD), ascorbate peroxidase (APX), catalase (CAT), glutathione peroxidase (GPx) and glutathione-S-transferase (GST) play important role in scavenging ROS in fungi under abiotic stress [42,43]. The study of gene expression involved in the alternation of antioxidant enzyme activities could provide insight into the molecular adaptation and toxicity of microbiology to heavy metal stresses. For instance, DEGs (including ggt, gpx, pepD, pepN) associated with glutathione metabolism were detected in Chryseobacterium sp. under Mt-Cd treatment [44]; Mn-SOD (sodA) accumulation in Escherichia coli P4 was exceeded by more than 2.6 folds after exposure to Cd²⁺ [23]; Markedly upregulation expression of antioxidant enzymes (SOD, POD, and CAT) in Pichia kudriavzevii, inhibited the Cd-induced ROS outburst and oxidative damage of proteins and membrane lipids [24]; The exposure of Chlamydomonas cells to the toxicants resulted in a progressive increase in mRNA transcript of Fsd 1, Msd 3 and Apx 1, up to a 2.5-fold increase compared with the levels detected in the control cells [45]; Twenty-four EpGST genes of Exophiala pisciphila (H93) were significantly up-regulated after exposing in 111.2 mg/L Cd²⁺ medium [25]. In our study, the expression of the enzymatic antioxidant genes was significantly regulated under three levels of Cd stress, including 5 Superoxide dismutase (SOD), 2 Catalase (CAT), 2 Ascorbate peroxidase (APX), 1 Peroxidase (POD), 5 Glutathione S-transferase (GST) and 2 Glutathione peroxidase (GPX) (Table S7). First of all, the expression of three CuSOD (DN51949_c4_g5, DN52158_c0_g2, DN53456_c0_g1) genes was highly upregulated and owed an ultra-high range value of log2FC (15.99–17.46). While the up-regulation rate of the MnSOD gene (DN52823_c1_g1) was on the low side. It indicated that CuSOD might have played a more important role in scavenging ROS in M. spongiola. Second, 2 CAT genes (DN54406_c0_g1, DN52979_c0_g4), 1 APX gene (DN47558_c0_g1) and two GPX gene (DN53754_c2_g4, DN51781_c0_g1) were only upregulated in 0.90 mg/L treatment and were not significantly expressed in other two groups, indicating that they might have only helped to scavenge ROS under 0.90 mg/L Cd stress. Third, one POD gene (DN54069_c0_g3), one APX gene(DN28932_c0_g1) and one GST (DN47015_c0_g1) were only upregulated in 1.50 mg/L treatment and were not significantly express under lower concentrations, suggesting that they might have only response to scavenge ROS more actively under higher concentration Cd stress. Fourth, the expression levels of four GST genes and two GSH-Px in M. spogiola under Cd stress were significantly up-regulated, except for GST1 (DN33532_c0_g1) which showed down-regulation, indicating that GST genes expression is considered to be closely helped to resistant of Cd stress in M. spongiola. Furthermore, only one GST gene (DN28953_c0_g2) was up-regulated in 0.15 mg/L Cd-treatment group, while most of the enzyme antioxidant genes were not differentially expressed at the concentration of 0.15 mg/L, suggesting that low concentration of Cd could not induce the expression of enzyme antioxidants. It was reported that MRP1 mainly transports structurally diverse amphipathic organic anions, most of which are conjugated with glutathione (GSH), glucuronide, or sulfate [46]. From this, we speculate that Morchella might dependent on a heavy metal transport system rather than an enzymatic antioxidant system to respond to Cd stress at low concentrations.

Thirdly, small molecular proteins such as metallothionein (MT), heat shock protein (HSP) and thioredoxin (TRX) usually act as non-enzymatic antioxidants in response to the presence of highly ROS in M. spongiola. The main roles of HSPs in fungi include folding, stabilization, transport and degradation of protein, which are involved in cell cycle progression, and they are also involved in the activation of many key signal transduction in fungi [26]. Hsp70 and Hsp12 play an important role in fungal morphogenesis under stress-related conditions [26]. The expression of the Hsp70 gene in Heliscus lugdunensis was time-dependent and dose-dependent, and it reached the maximum after 24 h of 80 mM Cd stress [47]. In our study, the expression level of Hsp70 was upregulated only in the 0.15 mg/L Cd-treatment group and Hsp12 was significantly upregulated in three Cd-treatment groups (Table S8). Combined with the morphological variation of M. spongiola mycelia under low concentration Cd stress (unpublished), it is speculated that Hsp70 and Hsp12 might responsible for maintaining cell morphogenesis and cell adhesion under Cd stress. Hsp90 plays a key role in promoting cell complexes that maintain genomic stability [48]. Hsp90 can trigger environmental adaptation under heat or cell wall stress by activating several key transcription factors and intracellular signal transduction proteins [49]. In M. spongliola, several Hsp90 genes were activated and expressed under cadmium stress, among which DN54532_c0_g10 was only up-regulated in 1.50 mg/L Cd-treatment group and the log₂FC value was as high as 16.29, indicating that Hsp90 plays an important role in the high concentration of Cd stress. Low molecular weight (LMW) Hsp is the most important protein produced by plants under heat stress by binding to partially folded or denatured substrate proteins to prevent irreversible aggregation, and there are few reports to prove that they are necessary for the function of normal cells [50]. In the transcript of M. spongiola under 0.90 mg/L Cd-treatment group, 4 LWM Hsp genes were significantly upregulated and no differential expression in the other two groups was occurred, indicating that LMW Hsp plays an important role in the stress response of M. spongiola under 90 mg/L Cd stress. Heavy metals can induce thioredoxin (Trx) messengers and change the redox state of cells, Trx-like genes can effectively increase the tolerance of bacterial to heavy metals [51–53]. Metallothionein (MT) binds to heavy metals through clusters of thiobonds [54]. The protein encoded by metallothionein MMT1 of Magnaporthe grisea has a high affinity for zinc and can be used as a powerful antioxidant [55]. In M. spongiola, 8 Trx were up-regulated and the log₂FC fluctuated within the range of 1.03–16.67. The expression level of five MT gene was ultra-higher (log₂FC range of 15.99–19.60) both in 0.90 mg/L and 1.50 mg/L treatments, indicating that Trx and MT might meritorious helped to increase tolerance ability in the M. spongiola antioxidation process under higher concentration Cd stress.

Forthly, Cd-responsive transcription factors have been found to play important roles in Cd stress response. Here, some transcription factor genes were significantly regulated in Morchella under Cd stress, including GATA, BTF3, MYB and SFC1 (Table S9). Basic transcription factor (BTF3) as known as a protein originally identified in HeLa cells that may be involved in the initiation of RNA polymerase II transcription [56]. Up-regulation of BTF3 could affect cell proliferation, apoptosis and cell cycle regulation of cancer cells [57]. It was also been demonstrated that BTF3 proteins were localized at the cytosol and nucleus of root cells and have positively regulated freezing tolerance [58]. MYB44 is controlled at multiple levels and strongly interconnected with MAPK signaling, its overexpression helps to increase the tolerance of Arabidopsis thaliana to osmotic stress and slightly improve the sensitivity to abscisic acid [59]. In our study, BTF3 was upregulated in 0.90 mg/L Cd-treatment group and the Myb44-like gene was upregulated in 1.5 mg/L Cd-treatment group, the log₂FC value of the above two genes was higher than 15.6. Therefore, it was speculated that all of the above transcription factors might play a critical role in M. spongiola to adapt to the different concentrations of Cd stress.

To summary, Figure 5 shows the Cd stress response regulatory networks in M. spongiola. On one hand, for the Cd environment, mycelia of M. spongiola may increase their Cd tolerance in multiple ways (Figure 5). Cd enters the cells mainly through the ion channels of HMT, ZRT, and CTR metal transporters, placing the cells under oxidative stress. This leads to the increased expression of genes encoding antioxidant enzymes including SOD, POD, CAT, and APX, as well as non-enzymatic antioxidants including TRX and HSP. Genes involved in glutathione metabolism are also up-regulated to further enhance ROS scavenging in the cell. Additionally, genes encoding metal chelating/binding proteins such as MT and HMT are up-regulated to sequester Cd and reduce its toxicity. Simultaneously, the expression of intracellular transcription factors including myb and btf is up-regulated, presumably affecting the activity of the related signal transduction pathways including the MAPK and cGMP pathways. Finally, Cd is excreted from the cells by upregulating the expression of membrane-bound metal transporters including MRP, CE, ILT and ABC. The combined effect of these mechanisms makes M. spongiola mycelia highly resistant to Cd stress. On the other hand, concentration-dependent response was persisted in RNA-seq analysis of the M. spongiola. The number of DEGs showed concentration-dependent increasing. DEGs in 0.90 mg/L and 1.50 mg/L Cd-treated group was 2.33 and 1.57 times higher than that in 0.15 mg/L group, respectively (Figure 1(A)). 0.15 mg/L could be defined as a tolerated concentration, where the number of differential genes including transporters, transcription factors, antioxidant enzymes and antioxidant non-enzymes was several times less than the other two concentrations (Tables S6–S9). 0.9 mg/L could be defined as cytotoxic concentration threshold, which mainly accompanied by differential expression of catabolism and cell cycle-related gene expression (Table 1). The results of this study also provide transcriptome level evidence of the differeneces how M. spongiola cope with Cd toxicity of different concentration.

Figure 5. Cd stress response regulatory networks in Morchella spongiola. ABC transporter: ATP binding cassette transporter; MRP: multidrug resistance protein; CTR: copper transporter; ZRT: zinc transporter; HMT: heavy metal transporter; CE: cation transport ATPase; ILT: Iron transport multicopper; SOD: superoxide dismutase; CAT: catalase; POD: peroxidase; APX: ascorbate peroxidase; GST: glutathione S-transferase; GPX: glutathione peroxidase; GSHD: glutathione dehydrogenase; HSP: heat shock protein; MT: metallothionein.

4. Conclusions

This study aimed to clarify the transcriptomic changes that occur in M. spongiola under Cd stress. The results obtained provide interesting insights into how M. spongiola cells perceive and respond to heavy metal stress. An RNA-seq analysis revealed 7444 genes exhibiting differential expression patterns in the mycelia of M. spongiola subjected to Cd stress. These genes were mainly involved in transcription-related processes, stress responses, antioxidant responses, and signal pathways. Based on these findings, a Cd stress regulatory network for M. spongiola was proposed. These results expand our understanding of the molecular mechanisms underpinning heavy metal tolerance in M. spongiola, and potentially in macrofungi more generally.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Natural Science Planning of Qinghai Province [2021-HZ-802] and Science and Technology Innovation Project of Qinghai University [2019-QH-16].