De novo transcriptome analysis reveals an unperturbed transcriptome under high cadmium conditions in the Cd-hypertolerant fern Athyrium yokoscense

Athyrium yokoscense shows strong tolerance to cadmium exposure, even at levels that are many times greater than the toxic levels in ordinary plants. To determine the mechanism of Cd tolerance in A. yokoscense , we grew these plants under high Cd conditions and observed the tissue-speciﬁc accumulation of Cd and generation of reactive oxygen species, which is one of the major physiological responses to Cd stress. Fuchsin staining indicated the existence of a casparian strip in A. yokoscense roots, which may participate in Cd hypertolerance in A. yokoscense . Moreover, we performed RNA-seq of RNA samples from A. yokoscense plants treated with or without Cd exposure and obtained comprehensive RNA sequences as well as the Cd-responsive expression patterns of individual genes. Through de novo transcriptome assembly and gene expression analysis, we found that A. yokoscense showed normal features with no signiﬁcant change in the expression levels of any transporter genes, even under high Cd exposure conditions. Our results demon-strate that A. yokoscense has an unusual mechanism that allows the invading Cd to partition into the distal roots, thus avoiding translocation of Cd into the xylem.


INTRODUCTION
Cadmium is very toxic to most plants when it is inadvertently taken up and translocated (Clemens, 2006). For example, Cd toxicity strongly affects nontolerant plants when their leaf concentration reaches 5-10 μg/g dry weight (DW) (White and Brown, 2010). Tobacco (Nicotiana tabacum) is sensitive to Cd stress, resulting in a significant reduction in plant weight along with severe inhibition of root elongation due to exposure to 100 μM Cd (Misra and Gedamu, 1989). In contrast, other plant species have acquired tolerance to high Cd concentrations and can grow after accumulating more than 100 μg/g DW Cd (Baker et al., 2000). One such species, Athyrium yokoscense, shows normal growth fea-tures under culture conditions containing 1,000 μM Cd and accumulates 1.5 to 3.5 mg/g DW Cd in the plant body over six weeks (Yoshihara et al., 2005). Athyrium yokoscense is a fern indigenous to eastern Asia including Japan (Fig. 1A) (Van et al., 2006). In addition to Cd tolerance, which reaches hundreds of times that of ordinary plants, A. yokoscense also displays tolerance to other heavy metal species (e.g., Zn, Cd, Pb and Cu) (Morishita and Boratynski, 1992). Thus, identifying the mechanism of Cd hypertolerance could advance our total recognition and understanding of metal tolerance and accumulation in plants.
Cd is often incorrectly absorbed into plant cells via specific transporters in the cell membrane that are ordinarily involved in the absorption of essential nutrient elements, such as Fe, Zn and Mn (Clemens, 2006). We have reported that the distal root tissues of A. yokoscense retained most of the absorbed Cd, but some could be translocated into the upper parts of the plant through the xylem (Yoshihara et al., 2014). This presumes one of the following two possibilities: (1) A. yokoscense has a specific transporter to separate essential nutrient elements from Cd, or (2) rhizome-like organs play a role in Cd tolerance in A. yokoscense. The rhizome is an underground organ in ferns that stores nutrients (Carlquist and Schneider, 2001). If the former is the case, little is known about the genomic information of fern species, including A. yokoscense. In fact, only the genome sequence of Selaginella moellendorffii, which belongs to the primitive fern clade Lycopodiophyta, has been determined thus far (Banks et al., 2011). If the latter is the case, it is important to know whether A. yokoscense contains a casparian strip among its underground parts. The casparian strip is mainly composed of lignin, fatty acid and suberin (Van Fleet, 1961), and restricts the influx of undesirable substances into the plant body (Steudle and Frensch, 1996). For example, in the roots, both water and nutrients are always transported across the casparian strip to vascular bundles (Steudle and Frensch, 1996). However, in ferns, there are few descriptions of the presence of the casparian strip (Lersten, 1997).
In this study, we observed the structure of A. yokoscense roots, particularly in relation to the casparian strip, and analyzed the effect of Cd and Cd-inducible stress on reactive oxygen species (ROS) in the roots. We also comprehensively examined A. yokoscense gene expression by de novo transcriptome analysis. Based on these results, we discuss fluctuations in gene expression and traits associated with hyperaccumulation of Cd in distal roots when the plants were treated with Cd. were prepared under in vitro conditions (Yoshihara et al., 2005). We used aseptically derived plants from the sporophyte stage that were redifferentiated from cultured cells (Fig. 1A). Athyrium yokoscense plantlets were grown in a magenta box containing 1/2 × Murashige and Skoog (MS) medium (1/2 diluted MS medium supplemented with 30 g/l sucrose, 0.3% gelangum, and vitamins) (Murashige and Skoog, 1962). The culture room was set at 25 °C, and plants were illuminated in a 14 h/10 h (L/D) cycle. Tobacco (N. tabacum L. cv. Petit Havana SR1) plantlets were used as a control, and were grown under the same conditions as the A. yokoscense plantlets.

Plant materials and growth conditions
Measurement of Cd content in plants Athyrium yokoscense was cultured for several months to the sporophyte stage, and N. tabacum was grown to a size similar to A. yokoscense (Fig. 1B, 1C). According to a previous report (Yoshihara et al., 2014) (Grabherr et al., 2013) with the parameters --seqType fq --SS_lib_type RF. Trinity-based contigs were further assembled using PCAP (Huang et al., 2003) with the parameters -y 10 and -t 80. To obtain a nonredundant transcripts dataset, PCAP-based contigs were grouped using CD-HIT-EST (Fu et al., 2012) to cluster sequences that have greater than 80% overlap with the longest sequence in a transcript cluster. Using the longest sequence in the transcript cluster, potential protein coding sequence regions were predicted using TransDecoder v.3.0.1 (https://transdecoder.github.io/) through sequence comparison with the UniProt database and Pfam-A database with a BLAST search (v2.3.0 + ) and the hmmscan program in the HMMER package (Johnson et al., 2010), respectively. Representative transcript sequences with coding potential were used for downstream analyses as the reference transcripts dataset. The gene-space completeness of the reference transcripts dataset was examined using BUSCO v3.0.2 (Simão et al., 2015) with -e 1e-05 -m tran and the embryophyta_odb9 database.
Transcriptome analysis The trimmed RNA-seq reads were mapped to the reference transcripts dataset using the BWA 0.7.17-r1188 program (Li and Durbin, 2010) with the parameters mem -M. The expression level of each transcript was quantified based on the reads per million mapped reads (RPM) computed from the read count data using the featureCounts program (Liao et al., 2014). Transcripts showing RPM ≥ 1 in all three samples were considered to be significantly expressed. Differentially expressed transcripts were identified using the DESeq2 program (Love et al., 2014) in R, with the thresholds as fold change ≥ 2 and adjusted q < 0.01.
Analysis of ROS and the casparian strip ROS in root sections were detected using CellROX Green reagent (Thermo Fisher Scientific) according to the manufacturer's protocol. The resultant roots were cut to 100 μm thickness sections using a tabletop hand microtome TH (Kenith, Osaka, Japan). They were observed using an upright microscope (ECLIPSE 80i, Nikon, Tokyo, Japan). Fluorescence was observed at 485 nm/520 nm (excitation/emission) with an ultrahigh pressure mercury lamp. Fluorescence intensity was measured using ImageJ image analysis software (W. Rasband, USA; http:// rsb.info.nih.gov/ij/). To observe the casparian strip, distal root sections were fixed in a formalin/acetic acid/alcohol mixture (3.7% formaldehyde, 5% acetic acid, 50% ethanol) for at least 24 h. The fixed samples were sequentially immersed in a series of ethanol mixtures at 50, 70, 85, 95 and 100%. Next, the root sections were embedded in Techno-vit 7100 (Kulzer, Hanau, Germany) and positioned. Sectioning was performed with a microtome PR-50 (Yamato Kohki Industrial, Saitama, Japan); the thickness of the paraffin ribbons was adjusted to 20 μm. The resultant ribbons were stained with 0.1% (w/v) New Fuchsin (Sigma Aldrich, St. Louis, MO, USA) and observed using an inverted microscope (PHASE CONTRAST-2 ELWD 0.3, Nikon).

RESULTS
Accumulation of Cd in roots and shoots and response to Cd stress To determine the response to Cd stress, we observed the morphology of A. yokoscense sporophyte-stage plants cultured under the low (0.1 μM) concentration of Cd condition and compared them with tobacco plants. There were no notable differences in appearance before and after Cd treatment in either species. In A. yokoscense, a large amount of Cd (4 μg/g DW) accumulated in the distal roots, whereas very low Cd accumulation (less than 1 μg/g DW) was observed in the proximal roots and shoots (Fig. 1D). In contrast, in N. tabacum, high levels of Cd accumulation (greater than 4 μg/g DW) were detected in both the distal and proximal roots (Fig. 1E). Thus, there was a significant difference in the Cd accumulation and its location in the tissues between A. yokoscense and tobacco.

De novo transcriptome assembly in A. yokoscense
We obtained RNA-seq reads using A. yokoscense sporophyte-stage plants cultured with or without 100 μM Cd. Using de novo assembly of these RNA-seq reads, we developed the first comprehensive nonredundant transcripts dataset of the Cd-hypertolerant fern. Via Illumina-based RNA-seq, we obtained 46 to 72 million paired-end reads for each sample (Table 1), and, in total,  (Table 2). We found that its N50 metrics were greater than 2.6 kb ( Table 2) and its transcript length distribution was largely greater than 2 kb ( Fig. 2A). We also assessed its assembly completeness using BUSCO (Simão et al., 2015), and found hits for 957 (66%) full-length and 70 (5%) fragmented BUSCO proteins (Fig. 2B)   searching against the NCBI nr database and found that 71.8% of the query sequences showed the highest hit to uncharacterized or hypothetical proteins (Fig. 2C), suggesting that the A. yokoscense genome encodes many functionally unknown genes. Moreover, we examined the distribution of A. yokoscense transcript orthologs in S. moellendorffii, M. polymorpha and A. thaliana using the OrthoMCL program (Chen et al., 2006), and found that 78% of A. yokoscense transcript orthologs are shared with at least one organism used to predict the orthologous groups (Fig. 2D).
Genes showing differential expression in culture with and without cadmium By comparing gene expression data between the samples with and without Cd exposure, we found an extremely small number of genes that were differentially expressed in response to Cd exposure (Fig. 3, Supplementary Table S1). When the analysis used the q-value as 0.00001, which was commonly employed, only four differentially expressed genes (DEGs), all without known functions, were detected in roots, and none was detected in shoots. When the q-value was set to 0.01, we found 43 and two DEGs in the roots and shoots, respectively (Table 3). Among these, all 43 transcripts increased in roots, while one gene increased and the other decreased in the shoots, under the Cd conditions (Table  3). These genes included a putative nitrate transporter (NPF family) gene (Athyo41034c000010) and three catalase genes (Athyo6890c000010, Athyo71933c000010 and Athyo38788000010), whose fold changes were calculated as 2.7-, 3.9-, 3.8-and 3.7-fold, respectively, in the roots (Supplementary Table S1). In shoots, their gene expression fold change values were not significantly different between treatment with and without Cd, suggesting that no alteration occurred in the shoots. No gene encoding another transporter was predicted among the significant DEGs other than the single NPF family transporter gene (Supplementary Table S1). We analyzed the expression levels of the NPF family gene and one of the catalase genes, because these expres- Fig.3. Ukai et al.  est DEG value of the three (Fig. 4, Supplementary Table  S1). The A. yokoscense transcripts dataset contained expressed gene sequences orthologous to genes reported to be involved in Cd transport (Supplementary Table  S2), such as NPF, AtIRT1, OsNramp5, OsHMA3, Hmt1, OsACA6 and AhCAX1 (Ortiz et al., 1992;Connolly et al., 2002;Miyadate et al., 2011;Ishikawa et al., 2012;Shukla et al., 2014;Ahmadi et al., 2018), whereas no ortholog for TaLCT1 (Clemens et al., 1998) was identified from this dataset using the BLASTP program. No significant changes were observed for the amounts of the following transcripts due to Cd stress: and another NPF gene (Athyo44698c000010) ( Table 4). Other orthologs were found, but they also showed no significant changes in expression levels (Supplementary Table S2).
Athyrium yokoscense generates few ROS in response to Cd stress and possesses a casparian strip Because exposure to heavy metals often causes strong stresses to plants via the generation of ROS (Andresen and Küpper, 2013), we analyzed ROS generation in A. yokoscense cells in the roots. ROS signals were detected in the central pillars and endothelial cells in the roots, but no marked difference was observed in the fluorescence intensity in cells with or without Cd exposure (Fig. 5A-5D).
Moreover, to understand the structural differentiation of A. yokoscense roots, sectioned roots were observed. Athyrium yokoscense roots have specific tissues corresponding to the epidermal cells and vascular bundles. We detected tissues stained with Fuchsin around the central pillars in root sections, in which large amounts of lignin were contained (Fig. 5E). Fuchsin is known to bind lignin and shows staining with red color (Barlow, 1969), and it can detect the presence of a casparian strip (Vermeer et al., 2014). These results indicate that there is a caspar-  sion levels increased in presence of Cd stress and showed greater than a 2.0-fold change in RPM. Quantitative RT-PCR results showed a slight increase in the amount of the NPF transcript, but no significant change in expression of the selected catalase gene, which showed the larg- ian strip in the endothelium.

DISCUSSION
In A. yokoscense, Cd mainly accumulated in the distal roots (Fig. 1D, 1E). In contrast, a small amount of Cd is transported to the proximal roots and shoots. We attempted to elucidate the molecular mechanism of tolerance to Cd in A. yokoscense cells. Since few data have been reported on A. yokoscense genes, we constructed an expressed sequence dataset through RNA-seq-based de novo transcriptome analysis of A. yokoscense cells with and without Cd exposure. We determined the RNA sequence details and expression patterns of individual genes in A. yokoscense. Our dataset was considered to have sufficient quality (Tables 1 and 2, Fig. 2). BUSCO analysis showed that our dataset contained contigs corresponding to complete cDNAs, and an Orthmcl analysis suggested that many predicted cDNAs in an A. yokoscense ortho group were shared with other organisms (Fig.  2). These results showed that this dataset is a representative database for cDNA sequences corresponding to mRNAs in common fern species.
Using these data, we analyzed the expression patterns that changed in response to Cd stress. Contrary to previously observed gene expression changes in normal plant species under Cd stress conditions, we found no significant changes in the A. yokoscense transcriptome, even under high Cd conditions (Fig. 3). It has been reported that expression of a large number of genes was altered in response to Cd stress in other plant cells. In the case of S. nigrum and S. torvum, the expression levels of 359 and 634 genes, respectively, were largely influenced when they were exposed to 50 μM Cd in the medium for 24 h (Xu et al., 2012). Similarly, a large number of DEGs have been reported in Oryza sativa under Cd stress: expression of 236 and 1,162 genes was changed by medium containing 10 μM and 100 μM Cd, respectively (He et al., 2015).
In contrast to these findings, we found that A. yokoscense showed no notable alterations in the expression levels of almost all genes in response to Cd stress (Fig. 3). We detected few genes showing differences in expression levels when the q-value was set to the commonly used value of 0.00001. Even when the q-value was set to a high value (0.01), we only found 43 and 2 DEGs in the A. yokoscense root and shoot transcriptome between plants with and without Cd stress (Table 4), indicating that the regulation of A. yokoscense genes was essentially unperturbed in response to Cd stress. Cd mainly moves into cells via transporters/channels, which show broad-range specificities for Fe, Ca and Zn ions (Clemens, 2006). To date, an IRT transporter belonging to the ZIP (ZRT-IRT-like proteins) family, an ortholog of the wheat LCT1 transporter, a cation (calcium) channel and Nramp have been reported to be involved in absorption of Cd into cells (Clemens et al., 1998;Connolly et al., 2002;White and Broadley, 2003;Ishikawa et al., 2012). We identified transcripts corresponding to orthologs of these transporter/channels in A. yokoscense, except for LCT1 (Table  3). These findings suggest that A. yokoscense has a similar ion transport system to other plants.
The Arabidopsis NPF, AtNPF1.8 (At4g21680.1), has been shown to be up-regulated by Cd exposure . In A. yokoscense, the expression level of an NPF ortholog (Athyo41034c000010) doubled in response to Cd stress (Fig. 4, and Supplementary Table S1). No significant changes were detected in the transcript levels of genes orthologous to other transporters/channels (Table  3). These observations imply that these genes are not involved in Cd absorption in A. yokoscense.
Catalase is an enzyme that catalyzes the reduction of H 2 O 2 to O 2 , and may be responsible for the removal of excess ROS during stress (Mittler, 2002). Oxidative stress is often induced by Cd exposure and causes severe toxicity (Andresen and Küpper, 2013). It has been reported that transgenic tobacco harboring a catalase gene, OsACA6, showed acquired tolerance to Cd stress (Shukla et al., 2014). Significant suppression of Cd-induced ROS is observed in an Arabidopsis helleri transformant overexpressing AhCAX1 (Ahmadi et al., 2018). In A. yokoscense roots, a slight elevation in the expression levels of catalase genes (Athyo6890c000010, Athyo71933c000010 and Athyo38788000010) was detected under the Cd stress condition by DEG analysis (Supplementary Table S1). However, no significant change was detected by qRT-PCR (Fig. 4B). DEGs of catalase genes showed less than 4-fold change (Supplementary  Table S1); they were detected by analysis using a high q-value. Therefore, it seemed that catalase gene expression was little influenced by the Cd stress. Balestri et al. (2014) reported that Petris vittata is tolerant to 60 μM Cd conditions, but its root structure is damaged when it is exposed to 100 μM Cd. In contrast, even when exposed to 100 μM Cd, A. yokoscense showed no visual symptoms (Yoshihara et al., 2005) and no significant ROS signals ( Fig. 5A-5D). We detected lower Cd concentrations in the proximal roots from A. yokoscense than in the distal roots (Fig. 1D, 1E). It has been reported that ROS are generated when Cd causes a deficiency in transition elements, such as iron, and reduces the role of glutathione working as an antioxidant (Andresen and Küpper, 2013). However, A. yokoscense generated little ROS under Cd exposure conditions (Fig.  5), whereas a large amount of ROS was generated in tobacco cells under Cd stress (Shukla et al., 2014). This suggests that Cd exposure causes low induction of ROS in A. yokoscense cells, implying that A. yokoscense has an unusual mechanism to avoid the translocation of Cd into the xylem.
Cross-sectional observation of the A. yokoscense roots showed the differentiation of tissues corresponding to epidermal cells and vascular bundles, as well as the existence of a casparian strip (Fig. 5E). A similar structure has been observed in the differentiated roots of P. vittata, another fern species, which has a casparian strip (Balestri et al., 2014). Interestingly, the casparian strip in A. yokoscense was stained by Fuchsin, suggesting that it possesses a lignified cell wall (Fig. 5E). In maize, it has been reported that strong lignification occurred in the pericycle and inner cortical cells under exposure to Cd stress (Lux et al., 2011). It is thus possible that lignification of the casparian strip is involved in tolerance to Cd stress in A. yokoscense.
Even under strong stress condition with a high concentration of Cd exposure, A. yokoscense maintains its normal features without a substantial change in transporter and channel gene expression levels. It has been shown that these genes are strongly influenced by Cd stress in other plant species (Xu et al., 2012;He et al., 2015). Nishizono et al. (1987) have reported that tolerance to heavy metals is largely due to the role of the root cell wall in A. yokoscense, as 70-90% of metal elements (Cd, Cu and Zn) are localized in the cell wall. Our results strongly support their observation and may explain why few genes showed marked expression changes due to Cd exposure.

AUTHOR CONTRIBUTIONS
Y. U. performed the major experimental work, analyzed the results and wrote the manuscript. H. S. conceived and planned the project and supervised the research. H. S., K. M. and T. Y. designed the experiments, discussed the results and drafted the manuscript. K. I. performed the bioinformatics analysis. M. K. and S. Y. performed the plant analysis under Cd stress. H. T., K. K., K. S. and F. G. contributed to writing the manuscript, discussing the results and revising the paper. All authors read the manuscript and agree with its contents.