Zostera japonica plants were collected from the intertidal seagrass beds (at low tide sites) from Huiquan Bay in Qingdao (E 120.21°, N 36.05°) in northern China in early April of 2015, and surface seawater samples from the sites were collected. Although the specific causes of the recent dramatic seagrass decline in Huiquan Bay are unknown, anthropogenic habitat deterioration (e.g., fishing) and environmental water contamination (e.g., wastewater discharge and port transportation) during the past decades of the twentieth century may be responsible. To ensure homogeneity, samples were vegetative propagated from a single genotype using stolon and rhizomes, and young new shoots were chosen for pot experiments.

Following a 24 h acclimation period, Z. japonica plants were planted in a hydroponic system in a growth chamber containing natural seawater from Huiquan Bay. The condition of the seawater was maintained as follows: salinity 32.3–33.6; temperature 19.3–20.4°C; dissolved oxygen (DO) concentration of 4.94–5.49 mg/L; oxidation-reduction potential (ORP) 44.1–50.6 mV; and a pH 7.84–8.01, natural light. At the start of the experiment, seawater used in the four sample replicates was supplemented with different concentrations (0 and 50 μM) of Cu, Pb, and Cd. Seagrass tissues were sampled seven days after the start of the treatments for RNA-seq (three replicates each) and qPCR (three replicates each) analyses; four replicates each for metal accumulation analyses.

Seagrass samples (four replicates) from all treatments at 0.5, 1, 2, 4, and 7 days were harvested and analyzed for metals' content. Details see the Methods section (Lin et al., 2016b).

Primers used to amplify differentially expressed RNA transcripts were designed for RT-qPCR analyses using Primer Express v. 3.0.1 software (Applied Biosystems, Foster City, CA, United States) (Table S4 in Supporting Information). RT-qPCR procedures outlined by Saminathan (Saminathan et al., 2015b) were used. The relative quantification based on three biological replications was normalized to the 18SrRNA reference gene and calculated using the method described by Livak et al Livak and Schmittgen (2001).

The heavy metal uptake in both above- and below-ground seagrass tissues increased with increasing concentrations of Cu, Cd, and Pb for all treatments (Figure 1). Initial days of Cd, and Pb showed no significant differences in metal uptake; and obvious differences only appeared only at Cu with increasing days; and detailed was higher than uptake in the roots. Regarding the Pb treatments, varied days (0.5–7 days) resulted in minimal variation between shoot and root. Furthermore, when exposed to Cd, the seagrass tissues accumulated more Cd as the treatment concentration and duration of exposure increased (Figure 1). Therefore, the typically treatment was Cu-contained, with higher uptake concentrations in both above- and below-ground seagrass tissues, and the uptake of the shoot from low to high concentrations was higher than that observed in the root (Figure 1).

Some gene regulation occurs in response poor conditions to which seagrasses are exposed. The quality and length distribution details associated with the raw sequence reads and clean reads were presented in Table S2 and Figure S1 (Supporting Information). Briefly, all 14.47 G of clean sequencing reads, derived from four libraries, were de novo assembled using the Trinity software. Gene expression levels were estimated using RSEM (Li and Dewey, 2011) for each sample. Clean data were initially mapped back onto the transcriptome, and the read count for each gene was then obtained from the mapping results. A full comparison of the DEGs for each treatment and a full list of the DEG levels and the associated FPKM density distribution for each treatment were shown in Figure S2 (Supporting Information).

As shown in the Venn diagram (Figure 2A), the overlapping regions (2001 individual genes) represent DEGs that were common in multiple treatments, and the non-overlapping regions represent DEGs that were unique to single treatments (Figure 2). In Cu, Pb, and Cd treatments, 1,862 (10.19%), 3,548 (19.42%), and 1,081 (5.92%) genes were differentially expressed (padj <0.05) (Figure 2A), respectively. Pb exhibited more DEGs than Cu and Cd, and 2001 (10.95%) genes were co-expressed in all three treatments. Furthermore, statistical analyses revealed that exposure to Cu, Pb, and Cd treatments caused the observed differential expression (padj < 0.05) of 6201 genes (3266 upregulated and 2935 downregulated), 7250 genes (4729 upregulated and 2521 downregulated), and 4815 genes (2388 upregulated and 2427 downregulated), respectively (Figure 2B). There were more upregulated than downregulated genes in all three metal treatments, which indicated that gene activation/upregulation was a common feature of the three metal treatments (Figure 2B and Table S3 in Supporting Information). Together, these results suggested that exposure to heavy metals at 50 μM strongly affected gene expression in Z. japonica.

To visualize the differences in the gene expression profiles for all three heavy metal treatments, the DEGs were analyzed using a one-way ANOVA, and the most significant genes (padj < 0.05) were subjected to clustering analyses. As shown in Figure S3 (Supporting Information), the gene expression profiles induced by Cu, Pb, and Cd were divided into two subclusters. The right subcluster consisted of samples exposed to control, Cd and Cu-spiked solutions, and the left subcluster included samples exposed to Pb treatments. These results indicated that characteristic gene expression profiles can be used to discriminate each of the metal treatment groups from the Pb and controls, and they can discriminate all three metal treatments from each other.

Twelve genes in random (three replicates) were examined and further analyzed using RT-qPCR (Table S4 in Supporting Information). Based on the Z. japonica gene expression levels as measured by RT-qPCR, compared to those measured using RNA-seq methods, the results indicate that all 12 genes exhibited different expression levels in three treatments (p < 0.001; Figure 3). A highly significant correlation (R² = 0.84) was observed between the RNA-seq and RT-qPCR data sets. Moreover, the trend associated with expression changes of these genes based on qPCR was the same as that detected by the RNA-seq method, which confirmed that the RNA-seq method generated reliable expression data.

An enrichment analysis was used to further identify the genetic markers and biological processes influenced by Cu, Pb, and Cd treatments. The results indicated that 1,862 Cu-specific, 3,548 Pb-specific, and 1,081 Cd-specific genes were categorized into 2,129, 2,141, and 2,033 biological processes, respectively.

GO enrichment analyses were used to further assess the DEG sets, and the results were summarized with GO slim terms (Figure 4). Regarding the molecular function category, structural constituent of ribosome associated with Cd and Pb exposure, kinase activity associated with Cu exposure were abundant terms, and the abundance of the general term “antioxidant activity” was ranked as follows: Pb > Cu > Cd. (Figure 4A–C) The “translation,” “protein metabolic process,” and “cellular protein metabolic process” terms were the top biological process terms associated with the Cu, Pb, and Cd treatments (Figure 4D). Moreover, “response to oxidative stress” term was also common; thus, the change in gene expression was likely due to stress caused by various abnormal factors associated with the metal stress conditions. The cellular component terms “cellular component,” “cell part” and “cell” were among the top abundant terms, which suggested the occurrence of changes in the activity of cells and the sequestration of heavy metal compounds in the vacuoles. In general, the abundance of these terms was ranked as follows: Pb >Cd >Cu. Furthermore, despite an increased number of genes/transcripts associated with Pb treatments, the results indicated that increased accumulation of metals under Pb conditions might result from a greater number of genes/transcripts in the transporter category as compared to those associated with Cd and Cu. These genes/transcripts are transporters, which play key roles in sequestering heavy metals into vacuoles without having detrimental effects on physiological homeostasis.

Our results indicated that the Cu concentration was higher in shoots than in roots, which possibly reflects the effective Cu ion exclusion mechanisms in the roots. The key mechanism of metal accumulation is the release of chelators from plant roots into the rhizosphere, which increases the solubility of metals in terrestrial plants. In aquatic plants, chelators might be released from both the shoots and roots, Z. japonica seemed to be a better hyperaccumulator of Cu, Pb, and Cd. Zostera.marina leaves,which can be a source of Cu, Pb, and Cd, and the leaves can cause Cu and Pb accumulation in sediments in seagrass ecosystems, but Cd were significantly lower in the surface sediment than in the leaves (Hosokawa et al., 2016). Using our cultivation conditions in seawater, high metal tolerance and hyperaccumulation in Z. japonica might be mediated by metal-pumps such as heavy-metal ATPase transporters (HMA) (Papoyan and Kochian, 2004), NAC TFs (N-terminal DNA binding domain transcription factors) (Liu et al., 2014), natural resistance-associated macrophage proteins (NRAMPs), and/or ABC transporters (Curie et al., 2000; Vert et al., 2002; Verret et al., 2004). ATPases as far as we known, which pump out metals using ATP to drive the reaction, and proton antiports, which use the proton gradient to pump metals across the cell membrane in the previous research (Nies, 2003).

The translocation factor (TF) was used to evaluate the ability of seagrass to transfer metals from the root to tissues above.

If the high enrichment plant transfer coefficient is >0.5 (Baker, 1981), then Cu, Pb, and Cd are all at high enrichment. This result likely occurs because aquatic plant leaves direct the absorption of heavy metals in the water. Results shows the increased Cu concentrations and transfer coefficients, which initially decreased before rising. Moreover, on the fourth day, TF peaked at the 50-μM concentration, and the value increased by 288% compared to the control concentration. Therefore, the order of heavy metal transfer sensitivity in Z. japonica was Cu > Cd > Pb (Table S1 in Supporting Information). The similar result also found in Zostera.marina leaves, which can be a source of Cu and Pb cycling in eelgrass ecosystems by accumulating after shedding and decomposing. But Cd were significantly lower in the surface sediment than in the leaves, the Cd concentrations that decreased during leaf decomposition increased in the surface sediment because of leaves that still contained higher Cd (Lee et al., 2004).

Our results indicated that the transcriptomes associated with the three metals dominated almost all of pathway categories, including Ribosome, Biosynthesis of amino acids, Carbon metabolism, Carbon fixation in photosynthetic organisms, Alanine, aspartate, and glutamate, Glyoxylate and dicarboxylate metabolism, etc. (Table 1). Metabolomic investigations have been conducted on terrestrial plant abiotic stress, and plant–pathogen interactions were recently reported (Kaplan et al., 2004; Allwood et al., 2008; Shulaev et al., 2008). Generally, the contributions of metabolic responses were found to respond more quickly to abiotic stress, and changes in transcript abundance from different biosynthetic pathways correlated with changing metabolite levels (Chechik et al., 2008; Ralser et al., 2009).

To test the hypothesis that the RNA-seq analysis is of Z. japonica is an effective tool and to explore the effects of the three metals on Z. japonica gene expression profiles, gene expression data were analyzed for co-expressed genes and genes with specific differences (i.e., Cu-specific, Pb-specific, and Cd-specific) under seagrass stress response was discussed and analyzed.

The genes that had the most significant expressed on Z. japonica were first identified. Therefore, genes that were co-expressed across different treatments may exhibit highest maximum likelihood values. Co-regulators associated with expression differed among all three metal treatments, as is shown in the Venn diagram (Figure 2A), and the overlapping regions (2001 individual genes) represent DEGs that were common in multiple treatments. The co-expressed genes were further clustered using DAVID functional annotation clustering at the highest stringency classification (Figure 5A). These process groups consisted of biological processes that have similar biological meaning because of shared genes. The results identified co-expressed genes that were associated with metal ion binding, response to chemical stimulus and wounding, photosynthesis, oxidoreductase activity, and peroxidase activity in molecular function, the transmission of messages between cells, and metabolic process, electron transport biological process (Figure 4D). The annotations of common genes in these blocks underwent NCBI non-redundant protein sequences (NR) database analyses and the most closely related species Zostera marina database (Table 2), and hits with maximum alignment lengths and sequence identity were identified as probable orthologs.

Furthermore, most of genes reflect the cellular response to stimulus, certain key gene regulators such as c6564_g2 response to oxidative stress, c53035_g1 and c63283_g1 response to chemical stimulus and wounding, c37417_g1 about nutrient reservoir activity, c21323_g2 in peroxidase activity and c57819_g2 which is involved in oxidoreductase activity were significantly highly expressed of following exposure to the three metals (padj < 0.05). More than that, genes such as c42083_g1 in the electron transport, c41999_g1 ABCs signal transduction, c6558_g1 and c2586_g1 in photosynthesis, light reaction and MT_S were significantly upregulated or downregulated following both three exposure (Figure 5A and Table 2). In higher plants, MT induction seems to be metal-specific. Plant MT genes were proved to be induced by many abiotic and biotic effectors, and play important roles in maintenance of homeostasis related to essential metal transition, detoxification of toxic metals, and protection against oxidative stress (Zimeri et al., 2005). In P. oceanica, MT mRNAs were increased by Cu and Cd, but not by Hg (Giordani et al., 2000), and this result is similar to what we observed in Z. japonica plants that were exposed Cu, Cd, and Pb groups. Genes including PSBS, PSAG, psaJ, psbA, psbD, etc. (Dattolo et al., 2014) are usually recognized as photosynthetic response genes, and they were detected in this study at significant expression levels. Heat shock proteins (HSPs) are often employed as a primary response to dampen heat stress effects, and play important roles in seagrass adaptability and resilience and (Gu et al., 2012; Franssen et al., 2014), our study also found HSP20 gene expression in all three metal exposure treatments. Therefore, it can be speculated that HSP genes are sensitive to some abiotic stresses. Because most of the non-model organisms bioinformatics resources for seagrasses are limited. There is a need for user-friendly and interactive web portals, whereby the user can perform annotation and mine for candidate genes, as well as genetic markers. Thus, sequencing of many more seagrass species like Z. japonica is necessary.

Changes in the regulation of a detoxification process including upregulation of genes coding detoxifying enzymes, such as cytochrome P450, UDP-glucose transferase, the ABC-transporter, as well as HSPs as a stress response gene to heavy metal (copper, Cu) was first shown in the seaweeds species Laminaria digitata and E. siliculosus (Ritter et al., 2008, 2014). In seagrasses, oxidative stress-protective genes have oftern been recognized and associated with light response (Dattolo et al., 2013), thermal stress (Franssen et al., 2011, 2014; Winters et al., 2011) and extreme environments (Lauritano et al., 2015). While the catalase gene has been reduced to a single copy in Z. marina (possibly due to reduced xylemcharacteristics of submerged plants), whilst all three types of superoxide dismutase remain (Olsen et al., 2016). Glutathione-related transcripts have been isolated and profiled in P. oceanica, Z. muelleri, and Z. marina (Massa et al., 2011; Lauritano et al., 2015; Olsen et al., 2016; Pernice et al., 2016). In protein, enzymes related to cellular stress and photosynthesis genes in which the studies of seagrass in response to differing light conditions also revealed the differential expression of genes for (Dattolo et al., 2013, 2014).

We also found a number of genes that were uniquely expressed in each metal (Datasheets S1–S4, Figure 5B). Genes for c21603_g1 about transposase, c35877_g1 in the Phenylpropanoid biosynthesis are up- regulated, c31516_g1 in integral to membrane, and c40473_g1 in the proteolysis and metabolic process were down-regulated and all of them are uniquely expressed in Cd treatment plants. Plants from the Pb treatment expressed genes are more than Cd and Cu treatment, such as chromatin-remodeling complex ATPase ISWI2(ISS)( c36097_g2) in nucleus position which is responsible for chromatin remodeling function, c32733_g1 which engaged in aminoacyl-tRNA biosynthesis, transmembrane transport genes (c34886_g2), c36004_g2 acted as transferase activity, c36201_g2 in protein kinase PCTAIRE and related kinases are all highly up-regulated.Similarly, genes associated with c38267_g3 in Cu ion transport process, c34401_g1 about Aryl sulfotransferase, c32383_g1 in Ubiquitin-conjugating enzyme and c33655_g1 in premnaspirodiene oxygenase were specifically expressed in the Cu treatment group.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.