H2S mitigates inhibition of proliferation of T. thermophila under Cd stress
Heavy metal pollutants caused toxic effects on ciliates, and the effect varied according to the bioavailable concentration and nature of the heavy metal [22]. An assay using the motile response of Tetrahymena pyriformis, gave a sensitivity better than 1 μM and a toxicity threshold to 7 μM for Cd [23]. Cd caused a dose-dependent decline in the viability of T. thermophila [24]. To understand the tolerance level of Cd for T. thermophila, the half maximal inhibitory concentration (IC50) value of Cd was determined, and it was calculated to be 30 µM for T. thermophila cells at 6 h culture (Fig. 1a). H2S alleviates Cd toxicity in plants. Exogenous H2S recovered Cd-induced growth inhibition in Brassica napus, Arabidopsis, and barley [25-27]. 70 µM NaHS (donor of H2S) largely stimulated proliferation of T. thermophila (Fig. 1b). Furthermore, the proliferation of T. thermophila under Cd and NaHS treatments was investigated to explore whether H2S physiologically alleviate Cd toxicity. An amount of 0.7×105 mL-1 cells was transferred to the SPP medium, and the number of cells was counted every 4 h. The proliferation inhibition of T. thermophila under 30 µM Cd was dramatically mitigated by 70 µM NaHS (Fig. 1c). The results showed exogenous NaHS play a protective role on Cd stress in T. thermophila.
Cd treatment promotes the production of endogenous H2S and cysteine in T. thermophila
In plant systems, endogenous H2S is generated through enzymatic pathways. Cysteine desulfhydrases regulates cysteine degradation into pyruvate, ammonia and H2S. However, O-acetylserine(thiol)lyase catalyzes the formation of cysteine using H2S and O-acetylserine. These physiological processes are interrelated under Cd stress [26]. Recently, we found cysteine is generated by reverse transsulfuration pathway involved CBS and CGL, and de novo pathway involved cysteine synthase (CS). At the same time, the CBS, CGL, and CS also catalyzed the H2S production in T. thermophila (in press) [28]. To explore whether endogenous H2S are involved in T. thermophila tolerance to Cd stress, formation of endogenous H2S were investigated. 10 to 30 μM Cd increased the H2S content and cysteine levels in a dose-dependent manner. When T. thermophila cells were treated with 40 or 50 μM Cd, both H2S and cysteine levels decreased due to stronger Cd toxicity (Additional file 1: Figure S1a-b). Exogenous cysteine treatment enhanced H2S level and maintained H2S at high level under Cd stress (Additional file 1: Figure S1c).
H2S alleviates lipid peroxidation and improves antioxidant capacity under Cd stress
Cd significantly inhibits the growth of microorganisms and plants. The treatment using 30 μM Cd also led to the stunted growth and nigrescence of T. thermophila after being exposed for 24 h. However, the toxic symptoms were drastically alleviated with NaHS supplement. The exogenous NaHS significantly inhibited Cd accumulation (decrease by 26%) in the T. thermophila cells (Fig. 2a). The H2S-mediated Cd accumulation was significantly decreased with hypotaurine (HT, a H2S scavenger that reverses the effect of H2S) treatment. NaHS application increased H2S level by 38% but did not affect it in H2S + HT group compared to untreated control (Fig. 2b). The results implied the regulatory role of H2S in Cd accumulation in T. thermophila.
NaHS increased the insoluble Cd fractions in salix leaves and roots [29]. Schizosaccharomyces pombe directly scavenge the free Cd2+ ions and the detoxification process occurs through the production of CdS nanoparticles [30]. In T. thermophila, spherical CdS nanoparticles in yellow colour with an average particle diameter of 186.9 ± 60.8 nm (Additional file 2: Figure S2a) were observed under Cd treatment, and the nanoparticles amount increased by adding NaHS (Additional file 2: Figure S2b and Additional file 3: Figure S3). However, UV-visible spectrum analysis showed that no CdS formation was found in vitro (Additional file 4: Figure S4). The formation of Ag nanoparticles from Ag ions was one of the defense mechanisms of T. thermophila against the toxic silver ions. Compared to AgNO3, Ag nanoparticles were remarkably less toxic. The Ag nanoparticles stored intracellularly in the food vacuoles of T. thermophila [31]. The results showed that the formation of CdS nanoparticles decrease Cd bioavailability and toxicity in T. thermophila.
Cd stress caused lipid peroxidation and induced malondialdehyde (MDA) production of T. thermophila cells in a dose-dependent manner (Additional file 5: Figure S5). Low concentrations of Cd had less effect on the MDA content in T. thermophila cells. However, 30 and 50 μM Cd lead to the MDA content of the cells increased by 91% and 395%, respectively. In comparison, the MDA content of the cells had no significant changes when the cells were treated with NaHS. However, the NaHS treatment markedly decreased the MDA content of T. thermophila cells under Cd stress (Fig. 2c).
It is well known that antioxidant defense system increases organism tolerance against metal-induced toxicity by upregulating the nonenzymatic antioxidants and different antioxidant enzymes. H2S increase GSH content and antioxidant enzymes activity in Arabidopsis [26]. Under Cd stress, GSH content increased by 51% and exogenous NaHS supplement further increased the GSH content in T. thermophila (Fig. 2d). But, the increase of GSH was reversed with HT supply. Furthermore, SOD activity increased by 92% and CAT activity increased by 29% under Cd stress. H2S supplement also enhanced SOD and CAT activities in T. thermophila cells (Fig. 2e, f). Combined application of NaHS and HT decreased the activities of SOD and CAT. The results indicated that H2S could alleviate Cd toxicity by improving the antioxidant capacity of T. thermophila cells.
Characterization of the sequenced Illumina libraries
Tetrahymena evolved various efficient detoxification pathways allowing the survival from Cd stress, such as overexpressing metal chelators and activating antioxidant signal pathways. The oxidative stress related mechanism of Ag nanoparticles was revealed at the transcriptional level. Some oxidative stress related genes were upregulated upon exposure to sub-lethal concentrations of Ag compounds, although intracellular ROS levels and SOD and CAT activities were not elevated in Tetrahymena [32]. To further explore the mechanisms of H2S alleviating Cd stress in T. thermophila, RNA-seq was employed to investigate the changes in genome-wide gene expression for four groups of cells: exposed to nutrients only (CK), with 70 μM NaHS (S), under 30 μM CdCl2 (C), and 70 μM NaHS+30 μM CdCl2 (CS), with three biological replicates. A total of 305.2 million pair-end reads with Q20 > 97% and Q30 > 93% were obtained from 12 libraries (Table 1). Among the short clean reads, more than 94% were mapped to the T. thermophila Functional Genomics Database (http://tfgd.ihb.ac.cn/). Approximately 95% of the reads from each library were perfectly matched to the reference genes, and more than 93% of the reads in the libraries were mapped to single locations. Half of these uniquely mapped reads in each library were mapped to the sense strand, whereas the other half was mapped to the antisense strand (Additional file 6: Table S1). Then, the mapped reads were further classified and annotated using TopHat [33]. The correlations between the three replicated samples were calculated on the basis of the normalized expression results (Additional file 7: Figure S6). The correlation coefficient between the three replicated samples was reasonable for the CK, S, and CS groups, but that between C1 and C2 or between C1 and C3 was lower than 70%. Thus, the C1 sample was abnegated.
Differentially expressed genes (DEGs) in response to NaHS treatment, Cd stress, and NaHS treatment under Cd stress
DEGs were hierarchically clustered to obtain a comprehensive view of the differential gene expression under NaHS treatment, Cd stress, and NaHS treatment with Cd stress (Fig. 3a). Under NaHS treatment, the expression level of 191 genes changed. Among them, 134 genes were upregulated and 57 genes were downregulated. Under Cd stress, the expression level of 9152 genes significantly changed, including 4658 upregulated genes and 4494 downregulated genes. A total of 1087 genes were upregulated and 272 genes were downregulated with NaHS treatment under Cd stress. The expression levels of most genes recovered under Cd stress with NaHS treatment. A total of 4122 genes were upregulated and 3738 genes were downregulated between Cd stress and NaHS treatment under Cd stress (Fig. 3b).
The 95 DEGs overlapped between CK vs. S (50.0%) and C vs. CS (1.2 %), indicating that H2S-responsive transcripts under normal condition were far fewer than those under Cd stress. The 806 DEGs in CK vs. C (8.8%) were common to CK vs. CS (59.3%), and the 6001 genes overlapped between CK vs. C (65.6%) and C vs. CS (76.3%) (Fig. 3c), suggesting most of the Cd-responsive transcripts were altered by H2S. The absolute value of the log2 ratio ranged from 1.00 to 14.70 in CK vs. C, and ranged from 1.00 to 12.21 in C vs. CS (Fig. 3d). The significantly upregulated genes under Cd stress (log2FC > 8) were considerably related to oxidoreductase, GPXs, GSTs, heat shock protein, and MTs (Additional file 8: Table S2). By systematic bioinformatics approach, the predicted T. thermophila Cd proteome included thioredoxins, heat shock proteins, GPXs, GSTs, and MT protein [34]. Compared with Cd stress, the unigenes significantly downregulated in the NaHS treatment under Cd stress (log2FC < -8) were mainly related to oxidoreductase, GPXs, GSTs, and heat shock protein (Additional file 9: Table S3). The results indicated that the redox system is sensitive for NaHS treatment and Cd stress in T. thermophila.
Gene ontology (GO) enrichment analysis of DEGs
To obtain the functional annotations of the DEGs for Cd stress and H2S treatment under Cd stress, GO category enrichment analysis was performed. For the comparison of CK vs. C, the 5740 DEGs were classified as 50 functional groups (Additional file 10: Figure S7). The functional groups were divided into three categories: biological process, molecular function, and cellular component. The biological process mainly comprises DEGs involved in metabolic process (2619, 45.63%), cellular process (2589, 45.10%), single-organism process (1406, 24.49%), biological regulation (777, 13.5%), and localization (706, 12.30%). In the category of cellular components, membrane (2711, 47.23%), membrane part (2517, 43.85%), cell (1798, 31.32%), cell part (1784, 31.08%), and organelle (1094, 19.06%) were the most represented groups. Among the molecular function category, the major groups were catalytic activity (3081, 53.68%), binding (2160, 37.63%), and transporter activity (416, 7.25%). For the comparison of C vs. CS, the 4883 DEGs were also classified as 50 functional groups. Between C vs. CS and CK vs. C, they had exactly identical classification patterns (Additional file 10: Figure S7).
Next, TopGO enrichment analysis was performed to obtain a detailed classification through false discovery rate (FDR) adjusted P-value of < 0.05 as the cutoff (Fig. 4). The distribution of enriched GO terms indicated that several DEGs were involved in oxidation–reduction process (GO:0055114), oxidoreductase activity (GO:0016491) and glutathione peroxidase activity (GO:0004602) in both CK vs. C and C vs. CS. Under Cd stress, 283 DEGs were included in oxidoreductase activity and 231 DEGs participated in oxidation-reduction process. Compared with Cd stress, 263 DEGs constituted the oxidoreductase activity with NaHS addition, and 182 DEGs involved in oxidation–reduction process. Furthermore, response to oxidative stress (GO:0006979) was enriched in CK vs. C. Cell redox homeostasis (GO:0045454) was enriched in C vs. CS. These data indicated that H2S responds to Cd stress mainly through the adjustment of the redox balance.
Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway enrichment analysis
The annotated T. thermophila transcripts were mapped to the KEGG pathways to investigate the genes involved in important metabolic pathways. Under Cd stress, the 1116 DEGs were mapped to the 252 KEGG pathways. For NaHS treatment under Cd stress compared with Cd stress, 966 DEGs were mapped to 247 KEGG pathways. The pathways considerably related to carbon metabolism, GSH metabolism, drug metabolism–cytochrome P450, and metabolism of xenobiotics by cytochrome P450 (Fig. 5a, b). Under Cd stress, 54 DEGs (4.84%) were distributed in the carbon metabolism, and 48 DEGs (4.97%) also enriched in this pathway with adding NaHS. The KEGG pathway of GSH metabolism includes primarily GPX and GST. 52 DEGs (4.66%) were enriched at the GSH metabolism under Cd stress, and 54 DEGs (5.59%) were also enriched with adding NaHS. Under Cd stress, 39 DEGs (3.49%) or 40 DEGs (3.58%) were distributed in drug metabolism–cytochrome P450 or metabolism of xenobiotics by cytochrome P450, and 40 DEGs (4.14%) or 42 DEGs (4.35%) were enriched in the two pathways with NaHS addition, respectively. Besides, under Cd stress, 67 DEGs (6.00%) were distributed in ABC transporters (Additional file 11: Figure S8), and 70 DEGs (7.25%) were enriched in this pathway with adding NaHS. The considerably enriched pathways in C vs. CS were similar to those of CK vs. C, indicating the important detox effects through H2S under Cd stress.
Regulation of cytochrome P450 and GSH metabolism under Cd stress and with NaHS treatment
Cytochrome P450 represents an important participant in regulatory networks of organism responses to Cd stress. The expression of genes encoding cytochrome P450 family proteins was strongly induced by Cd in rice [35]. In the wolf spider Pardosa pseudoannulata, cytochrome P450 genes were found to respond to Cd stress [36]. A total of 44 putative functional cytochrome P450 genes were identified and classified into 13 families and 21 sub-families according to standard nomenclature in T. thermophila. Cd induced the expression of 6 cytochrome P450 genes and inhibited the expression of 9 cytochrome P450 genes (|Log2FC|>1). NaHS addition reversed the expression of 11 of them (Table 2). Under Cd stress, metabolism of xenobiotics by cytochrome P450 was significantly enriched by KEGG pathway analysis on DEGs, of which 32 DEGs were markedly upregulated (Additional file 12: Figure S9), including 30 GSTs, and 29 of them were markedly downregulated when NaHS was added.
63 cytosolic GSTs have been identified in T. thermophila. Under Cd stress, GSH metabolism was significantly enriched by KEGG pathway analysis on upregulated genes (Fig. 6). The upregulated DEGs included 31 GSTs and 8 GPXs. Compared with Cd stress, GSH metabolism was also significantly enriched on downregulated genes with NaHS addition, including 30 GSTs and 8 GPXs. Specifically, the 8 upregulated GPXs and 29 GSTs induced by Cd stress returned to normal level with adding NaHS. GSTs detoxify xenobiotic compounds by linking the -SH group of antioxidant GSH covalently to a substrate [37]. Given the role of GSH in the detoxification of heavy metals, regulation of GSH, GPXs, and GSTs by H2S play an important role in ameliorating Cd stress.
Regulation of ABC transporters under Cd stress and with NaHS treatment
GSH marked xenobiotic molecules are excreted from cells via phase III proteins such as ABC transporter enzymes. Based on the KEGG pathway analysis, the ABC transporters related to detoxification of Cd in T. thermophila was screened and listed in Fig. 7. Under Cd stress, 49 DEGs involved in the enrichment of ABC transporters were markedly upregulated, while 18 DEGs were downregulated (Fig. 7). The ABC transporter family of T. thermophila was classified into 8 distinct groups [18]. The ABCA subfamily comprises of a total of 32 members, which included 9 upregulated genes and 3 downregulated genes under Cd stress. The ABCB subfamily consists of 26 transporters, of which 5 genes were upregulated and 5 genes were downregulated. The ABCC subfamily comprises 60 transporters and represents the largest family of ABC transporters in T. thermophila. 22 genes were induced by Cd and 6 genes were downregulated. The ABCG subfamily includes 39 transporters, and 9 genes were upregulated under Cd stress and 3 genes were inhibited. ABCG19 (TTHERM_00034920) was induced up to 357 folds under Cd stress (Additional file 3: Table S2), which suggested it could play an important role in detoxifying Cd toxicity. H2S has been demonstrated for direct regulating ABC transporters to induce stomatal movement in plants. Exogenous H2S induces stomatal closure and this effect is impaired by the ABC transporter inhibitor glibenclamide in guard cells [38]. Compared with Cd stress, 42 DEGs were downregulated and 28 DEGs were upregulated for ABC transporters with NaHS addition. Most of the expression of ABC transporters under Cd stress was restored and some of them even hyper activated by exogenous NaHS addition, which is in relation to H2S lowering the cellular concentration of Cd (Fig. 2a) and promoting CdS nanoparticles formation (Additional file 2: Figure S2b and Additional file 3: Figure S3).
Quantitative real-time polymerase chain reaction (qRT-PCR) validation of differentially expressed transcripts
To validate the accuracy and reliability of transcriptome sequencing, six DEGs involved in chelating heavy metals, redox reaction, and stress response were evaluated through qRT-PCR. NaHS slightly affected the expression of the genes under non-stressed condition. Under Cd stress, the expression levels of the genes were significantly upregulated. The log2FC of the CdMT genes were ranked as MTT5≈MTT3>MTT1 under Cd stress (Additional file 3: Table S2). The expression levels of four genes decreased with NaHS addition, including SSA6, OXR1, GLR2, and MTT3. On the contrary, the expression level of MTT1 further increased and MTT5 had no obvious change with NaHS treatment under Cd stress (Fig. 8). The three CdMT isoform genes present differential expression patterns between Cd stress and NaHS treatment under Cd stress, indicating their functional diversification. The correlation between qRT-PCR and RNA-seq was measured by scatter plotting log2-fold changes, which showed a positive correlation coefficient in both techniques (Pearson coefficient R2 = 0.95), thereby indicating the reliability of the sequencing data.