Characterization of rhizome transcriptome and identification of a rhizomatous ER body in the clonal plant Cardamine leucantha

The rhizome is a plant organ that develops from a shoot apical meristem but penetrates into belowground environments. To characterize the gene expression profile of rhizomes, we compared the rhizome transcriptome with those of the leaves, shoots and roots of a rhizomatous Brassicaceae plant, Cardamine leucantha. Overall, rhizome transcriptomes were characterized by the absence of genes that show rhizome-specific expression and expression profiles intermediate between those of shoots and roots. Our results suggest that both endogenous developmental factors and external environmental factors are important for controlling the rhizome transcriptome. Genes that showed relatively high expression in the rhizome compared to shoots and roots included those related to belowground defense, control of reactive oxygen species and cell elongation under dark conditions. A comparison of transcriptomes further allowed us to identify the presence of an ER body, a defense-related belowground organelle, in epidermal cells of the C. leucantha rhizome, which is the first report of ER bodies in rhizome tissue.

www.nature.com/scientificreports/ Thus, the development of rhizomes and the controls that regulate this process have been extensively studied in rhizomatous crops and wild plants 7,9,10 . A comparison of transcriptomes between organs/tissues at the different developmental stages resulted in the detection of contrasting gene expression patterns [10][11][12] . Furthermore, unique genes expressed in particular organs/tissues have been identified 7,13 . Transcriptomic studies of rhizomes have been conducted in species such as Sorghum sp., Phragmites australis and Oryza longistaminata 7,9,10,14 . These studies have identified genes specifying rhizomes by comparing gene expression between rhizome tissues, e.g. tips and internodes 7,9,10,15 , between rhizomes and above-ground shoots or leaves [16][17][18] and between rhizomes, roots and above-ground organs 16,19,20 .
A comparison of the transcriptomes of shoots and rhizomes derived from the same meristem but exposed to different environments under natural conditions provides concise information about the developmentally controlled robustness and plastic responses of plant shoot transcriptomes to the environment. Although developmental identity is considered to be the strongest determinant of transcriptomes, gene expression should also be largely determined by the environment surrounding an organ 21,22 . Differences in gene expression between rhizomes and shoots should therefore partly reflect the difference in physical and biological environments above and below ground. Furthermore, comparative transcriptomic analysis is expected to identify whether a rhizome has a specific system that is known to be either shoot or root specific. For example, the ER (endoplasmic reticulum) body is a log-shaped organelle constitutively present in the root and hypocotyl of Arabidopsis thaliana and related plants (predominantly in the Brassicaceae family) and contains a large amount of proteins with myrosinase (β-thioglucosidase) activity 23,24 , which is crucial for the production of key defense substances of Brassicaceae, i.e. glucosinolates in response to attack by herbivores and microbial pathogens. Thus, it has been proposed that ER bodies constitute the defense machineries of the roots, hypocotyls and cotyledons of A. thaliana 25,26 . Although the Brassicaceae family contains many rhizomatous species, it remains unclear whether ER bodies are present in rhizomes, one of the major belowground structures of plants.
In this study, we aim to characterize the rhizome transcriptome of Cardamine leucantha (Brassicaceae) under natural conditions. The species is a clonal herbaceous plant that grows on deciduous forest floors and along forest margins. An aerial shoot elongates aboveground in spring from the tip of a rhizome that is produced in the previous year. The sequential expansion of leaves from aboveground shoots is followed by the production of flowers and the growth of the aboveground shoots stops with flowering and fruiting ( Fig. 1a). At the same time, C. leucantha produces belowground, creeping, stoloniferous rhizomes at the base of the growing shoot (Fig. 1a). According to the morphological definition 6, the C. leucantha rhizome is classified as a secondary rhizome, originating from the lateral meristem of the main axis that forms the aboveground shoot. The developmental patterns of C. leucantha allowed us to simultaneously harvest four distinct tissues/organs, i.e. aboveground shoot, leaf, rhizome (belowground, derived from the shoot meristem) and root (belowground, derived from the root meristem), growing under identical climate conditions in the spring. In particular, we aimed to characterize the rhizome by comparing its transcriptome with other representative tissues. We specifically asked the following questions: (1) How similar is the transcriptome of the rhizome to that of the aboveground shoot and does their similarity reflect a shared developmental origin? (2) How similar is the transcriptome of the rhizome to that of the root and does their similarity reflect their shared belowground environment? (3) Are there any genes that are specifically expressed only in rhizomes?

Results and discussion
Gross morphological and anatomical observations of rhizomes. We conducted observations and sampling in a natural population of C. leucantha at a site in Rikubetsu, Hokkaido, Japan (43°27′ N, 143°46′ E, 251 m a.s.l.), located within a cool-temperate deciduous forest along the Toshibetsu River. The forest was dominated by Salix sachalinensis and contained Fraxinus mandshurica, Quercus crispula and Ulmus davidiana as common tree species. At this site, C. leucantha ramets elongated to form 30-to 60-cm upright stems and produced inflorescences with insect-pollinated white flowers in June. At the same time, one or more rhizomes started to elongate from the lateral shoot meristem at the basal part of the flowering shoot (Fig. 1b). A stoloniferous rhizome grew horizontally by elongating internodes (Fig. 1c), which then fully expanded by the end of the summer (23 cm on average, up to 65 cm in our observation). During the winter, the rhizome remained beneath the ground surface, when leaves and inflorescences had already formed at its tip. These shoot tissues appeared aboveground in the next growing season (Fig. 1d) and new rhizomes were observed at the base of the new shoot (Fig. 1e). We observed median longitudinal sections of tips of young rhizomes with toluidine blue stain under a microscope ( Fig. 1f-g). The rhizome apex showed a typical shoot apical meristem structure with tunica/corpus organization (Fig. 1g), consistent with its developmental origin as a shoot apical meristem.
Sequencing summary. We then collected RNA samples from C. leucantha plants at the study site in Rikubetsu to obtain the transcriptomic landscape of the C. leucantha tissues under natural biotic and abiotic conditions. We first collected RNA samples from four tissue types, i.e. a rhizome apex, a flowering shoot apex, a root apex and a leaf, at five time points from May 2011 to February 2012 for de novo assembly of the reference sequence with a high gene coverage. The cDNA sequencing data of these samples, obtained using the 454 Titanium platform (Roche, Basel, Switzerland), contained 1.5 M reads with an average length of 432 bp. In the final assembly using Newbler (454 Life Sciences; version 2.6), 27,834 isotigs (transcripts) were obtained with an average length of 1,386 bp and an N50 of 1,618 bp. These isotigs were then queried using the basal local alignment search tool (BLAST) Blastx (version 2.2.26) against Arabidopsis thaliana (TAIR10) protein data. 26,035 sequences (93.5%) successfully matched the database sequences with an e-value ≤ 1e−10. All 27,834 isotigs were used as the reference genes for the subsequent transcriptomic analysis. For transcriptomic resequencing, four issue types described above were collected from two plants in the flowering stage at noon on 31 May 2012 www.nature.com/scientificreports/ successfully mapped except for roots A and B (42% and 51%, respectively). The lower mapping rates in roots were partly due to the contamination of other organisms such as fungi, bacteria and virus, probably because the root samples were collected form natural forest soils.
Overall characterization of the rhizome transcriptome. To compare the overall gene expression patterns of rhizomes, shoots, roots and leaves, PCA (principal component analysis) was conducted. The first and second axes explained 35.2% and 26.7% of the total variance, respectively (Fig. 2a). The patterns of plants A and B were consistent in which the rhizome, shoot and root samples were well separated. Shoot and leaf samples were closely located according to both the first and second axes. Along the first axis, rhizomes were located between shoot and root samples (Fig. 2a). Thus, the overall gene expression patterns of the rhizomes shared characteristics with those of both shoots and roots to some extent. Along the second axis, the scores of rhizomes deviated from those of the other tissues. This indicated that the expression of genes explained by the second axis was unique to the rhizomes. The transcriptomes of the four tissues of C. leucantha were compared with microarray data from representative tissues of A. thaliana, i.e. hypocotyl, roots (7 days and 17 days old), vegetative shoot apex, rosette leaf, flower stage and pollen (AtGenExpress). The similarities of transcriptomes were evaluated by Spearman's rank correlation coefficients (Fig. 2b). For both plants A and B, the expression patterns of the rhizome of C. leucantha showed a relatively high correlation with those of the roots (7 and 17 days old) and hypocotyl and shoot apex of A. thaliana (Fig. 2b). The results support the idea that the rhizome shares characteristics with the shoot and root. The gene expression patterns of the shoot, root and leaf samples of C. leucantha were highly correlated with the corresponding tissues of A. thaliana. The shoot transcriptome of C. leucantha showed the highest similarity with that of A. thaliana flower likely because the shoot samples of C. leucantha in spring contained young flowering buds (Fig. 2b).
Transcript clustering based on tissue-dependent expression patterns. K-mean clustering resulted in classification of all expressed genes into 16 clusters based on the gene expression patterns across the four tissues from plants A and B (Fig. 3). The number of transcripts included in each cluster ranged from 42 (cluster 1) to 1,696 (cluster 5). Significantly enriched gene ontologies (GOs) within clusters ranged from seven in clusters 1 and 13 to 380 in cluster 5 (P < 0.05, Supplementary Table S1).
Remarkably, none of the clusters exhibited rhizome-specific expression patterns, whereas several clusters of genes were expressed in shoot-and root-specific manners (Fig. 3). Transcripts that showed relatively high www.nature.com/scientificreports/ expression in rhizomes were grouped in clusters 1-9 which were also highly expressed in the roots (clusters 1-3) or in the aerial tissues (cluster 4, Fig. 3). Cluster 1 was significantly enriched in cell wall-related GOs (GO:0015928, fucosidase activity; GO:0005199, structural constituent of cell wall) and symbiosis-related GOs (GO:0009610, response to symbiotic fungus; GO:0009608, response to symbiont, Supplementary Table S1). Clusters 5-9 were universally expressed across the four tissues. These observations were consistent with the intermediate locations of the rhizome transcriptome relative to the root and aboveground transcriptomes in the PCA plot. In contrast, the transcriptomic profiles of the other three tissues were characterized by the clusters specific to each tissue. Clusters 10 and 11, 12 and 13, and 14 and 15 were specific to leaf, root and shoot tissues, respectively and cluster 16, which was highly enriched by cuticle-wax-related GOs, represented transcripts that were specific to aboveground tissues (Fig. 3).
Representative transcripts of the rhizome as a belowground structure. Clusters 1-3 represented transcripts whose expressions were higher in the rhizome and root than in the leaf and shoot. In particular, genes in cluster 1 were exclusively expressed in rhizomes and roots and were nearly absent from the aboveground transcriptomes ( Fig. 3). Representative transcripts in cluster 1 contained those annotated to A. thaliana genes related to belowground defense and the control of reactive oxygen species (ROS), which suggested a functional aspect of the rhizome as a belowground organ. The top seven transcripts with the highest average rhizome expressions in cluster 1 (Table 1)  . PR-4 is involved in JA (Jasmonic Acid)-dependent defenses and is upregulated by treatment with necrotrophic fungi 27 and rhizobacteria 28 . A homolog of LSH10 in potato (Solanum tuberosum L.) has been reported to show very high transcript levels in stolons (rhizomes) and young tubers 29 . AtPrx37 encodes a peroxidase superfamily protein that has been reported to be expressed in the roots of A. thaliana, as well as in the basal parts of flowering stalks and mature rosette leaves 30 . Furthermore, AtPrx37 was reported to be up-regulated in A. thaliana roots in response to increasing ROS concentration under nitrogen, phosphorus, and potassium deficiency and is responsible for mineral uptake 31 . AtPrx39 was also reported to be involved in the control of the balance between distinct classes of ROS in the roots of A. thaliana, thereby regulating root meristem homeostasis 32 . MSRB9 has been reported to be prevalently expressed in A. thaliana roots and is involved in tolerance to the accumulation of ROS 33 . Overall, these findings suggest that the rhizome defense machinery possesses a belowground characteristic. Cluster 16, which contained genes that were specifically expressed in aboveground tissues but not in roots and rhizomes, also characterized the transcriptome of the rhizome as a belowground organ (Fig. 3). For instance, CER1 in this cluster encodes an enzyme that converts leaf/stem wax C30 aldehydes to C29 alkanes 34,35 suggesting that the composition of cuticular wax in rhizomes may be different from that of stems and leaves.
Representative transcripts of the rhizome as a shoot-derived structure. In contrast to the aforementioned clusters, transcripts in cluster 4 were highly expressed in rhizomes as well as in shoots and www.nature.com/scientificreports/ leaves, but only moderately in roots ( Fig. 3 and Table 2). The top seven transcripts with the highest average rhizome expressions among annotated genes in cluster 4 ( Table 2) were Lipid Transfer Protein 2 (LTP2, AT2G38530), β-glucosidase 18 (BGLU18, AT1G52400), a cell wall protein (AT2G10940), Epithiospecifier Protein (ESP, AT1G54040), Glutathione Peroxidase 3 (GPX 3, AT2G43350), Cell Wall-Plasma Membrane Linker Protein (CWLP, AT3G22120) and GAST1 Protein homolog 4 (GASA4, AT5G15230). The LTP2 transcript in A. thaliana was reported to be highly accumulated in the epidermal cells of the hypocotyl and cotyledons in dark-grown seedlings 36 . In A. thaliana, ESP was found to be consistently present in the epidermal cells of all aerial parts 37 and encodes a myrosinase cofactor, which is necessary to drive the myrosinase-catalyzed hydrolysis of glucosinolates and prompts terminal production of nitriles and epithionitriles in Brassica and Arabidopsis 38 . Collectively, high levels of transcripts homologous to these genes in C. leucantha rhizomes likely represent a shoot-derived tissue, particularly one growing under dark conditions.
Transcripts with high expression in the rhizome relative to other tissues. Because there was no k-mean cluster that showed strong rhizome-specificity in its expression, we compared the expression between four tissues for each transcript. We found that 394 transcripts annotated to 172 A. thaliana genes showed the maximum expression level in rhizomes, with twofold higher expression compared to the tissue with the second highest expression level (Supplementary Table S2). The top three genes in the difference of expression compared with the tissue with the second highest expression level showed 27-, 14-and 12-fold differences in the FPKM value and 23 transcripts showed more than fivefold differences (listed as the Log 2 FPKM difference in Table 3,  Supplementary Table S2). The top one was annotated to a gene (AT4G22485) encoding a lipid-transfer protein whose function has not been addressed ( Table 3). The transcript with the second highest expression level was annotated to a gene coding a cell-wall loosing protein, Expansin 3 (EXPA3, AT2G37640), that promotes cell  www.nature.com/scientificreports/ expansion in the roots of A. thaliana 39,40 . Furthermore, two additional transcripts within the top 10 were annotated to A. thaliana TUB1 (β-tubulin; AT1G75780) and ARR3 (type-A response regulator; AT1G59940), both of which are potentially related to root/hypocotyl elongation. TUB1 expression has been reported to be high in etiolated seedlings of A. thaliana and is suppressed by light 41 . ARR3 has been reported to be constitutively expressed in the vascular tissue of both shoots and roots and is induced by cytokinin in root tissues 42 . These observations are consistent with the extensive elongation of C. leucantha rhizomes under dark belowground conditions.

Enrichment of ER body-related genes and identification of ER bodies in the rhizome.
We noted that the transcripts annotated to the key genes for ER-body formation were highly expressed in the rhizome tissue of C. leucantha (Fig. 4a). A gene homologous to BGLU23/PYK10, encoding the major component of the ER bodies, was highly expressed in the rhizomes and roots of C. leucantha but was absent in the leaves and shoots ( Fig. 4a and Table 1). The homologs of genes encoding GDSL lipase-like protein 23 (GLL23, AT1G54010) and PYK10-binding protein1 (PBP1, AT3G16420) also showed high expression in both rhizomes and roots. In A. thaliana, these proteins have been reported to be members of the PYK10 complex 43,44 , which is a huge protein aggregate that facilitates the enzymatic activity of PYK10 45 . For transcripts homologous to the gene essential for ER body formation, NAI2, and its homolog TSA1 (AT3G15950 and AT1G52410) 46,47 , the gene encoding the major component of leaf-type ER bodies 48, were also highly expressed in rhizomes of C. leucantha ( Fig. 4a and Table 2). Furthermore, transcripts related to the biosynthesis of indole glucosinolate (CYP79B2, AT4G39950; CYP83B1, AT4G31500; TSA1, AT3G54640), which was proposed to be the in planta substrate of PYK10 24, were also enriched in the rhizome (Clusters 3 and 4,Supplementary Table S2). Enrichment of expressions of ER body-related genes in the rhizome prompted us to test the presence of ER bodies in C. leucantha rhizomes. We visualized the ER in the epidermal cells of rhizomes and leaves by transiently expressing green fluorescent protein (GFP) fused with a signal peptide and an ER-retention signal (SP-GFP-HDEL). In the rhizomes, rod-shaped structures were detected in addition to the typical ER network (Fig. 4b). These were absent in the leaf epidermal cells (Fig. 4c). These structures were similar to the ER bodies in A. thaliana in both shape and size, suggesting that they correspond to the ER bodies of C. leucantha. This observation was further corroborated by the electron micrographs of rhizomes that delineated the presence of a spindle-shaped compartment with ribosomes on its cytosolic surface, which resembled the ER bodies in A. thaliana. Interestingly, we found amyloplasts in rhizome cells, i.e. organs containing starch grains in a plastid, Table 3. Ten annotated transcripts of Cardamine leucantha that showed the highest expression in the rhizome relative to the tissue with the second highest expression. Cluster number in k-mean analyses (Cl. no.), transcript ID in C. leucantha (isotig number denoted by Newbler), gene expression (Log 2 FPKM) in the leaf, shoot, rhizome and root of plants A (A) and B (B), tissue with the 2nd highest FPKM value, FPKM difference between the rhizome and the tissue with the 2nd-highest expression (Log 2 FPKM), and AGI code and short description of corresponding annotated genes in the Arabidopsis thaliana database (TAIR10 protein sequences) are listed. www.nature.com/scientificreports/ but not chloroplasts (Fig. 4d), suggesting that the rhizome of C. leucantha serves as the energy storage organ [49][50][51] .

Cardamine leucantha transcripts
Overall, these results clearly demonstrate that the rhizome of C. leucantha develops ER bodies similar to the roots of the other Brassicaceae plants. Notably, in C. leucantha rhizomes, the ER bodies appeared to contain proteins homologous to PYK10 or to BGLU18, which are the major components of constitutive and facultative ER bodies, respectively. This might suggest that the C. leucantha rhizome ER bodies might represent both types of ER bodies, which is consistent with our comparative transcriptomic analysis that revealed the characteristic of rhizomes to be intermediate between the above-and belowground tissues.

conclusions
A comparison of the transcriptomes of rhizomes, shoots, roots and leaves in the rhizomatous plant C. leucantha revealed that the rhizome has a transcriptome characterized as an intermediate between the shoot and the root. Previous studies investigating rhizome transcriptomes have reported a clear difference in gene expression pattern between the rhizome and aerial shoots/leaves in bamboo 15, Oryza 16, Sorghum 7,52 and Atractylodes 18, however, its relationship to the transcriptome of belowground organs such as roots remained unclear. We included root samples in our analyses and showed that the rhizome and the root strongly share transcriptomic patterns likely due to the shared belowground environments. The overall characteristic of the rhizome transcriptome is intermediate between those of the aboveground shoot and the belowground roots, suggesting that both endogenous developmental and external environmental factors are important in the regulation of the rhizome transcriptome. Comparisons between rhizomes, roots, leaves and shoots have been reported for clonal grass species www.nature.com/scientificreports/ such as Oryza 16,53 and Miscanthus 19,20 . The reported lists of the genes that were highly expressed in the rhizomes in these studies did not correspond to that of this study, at least for the top twenty genes, probably reflecting distant phylogenetic relationships between monocots and dicots. The comparison of transcriptomes further allowed us to identify the presence of ER bodies, defense-related belowground organelles, in epidermal cells of C. leucantha rhizomes, which is the first report of an ER body in rhizome tissue. Our study suggested that the surrounding environment largely influenced organ identity at the gene expression level, but further studies are required to determine how much of the rhizome-specific expression patterns are constitutive and/or responsive to surrounding environments. De novo assembly samples were processed by the 454 Titanium platform to obtain cDNA sequences of C. leucantha. These de novo assembled data were used to construct the references onto which the RNA-Seq data were mapped. Total RNA was collected from four tissues in the different seasons described above using RNeasy Mini Kit (QIAGEN) before being purified twice by Oligo (dT) 25 Dynabeads. All samples were pooled as a 290ng mRNA sample to synthesize cDNA. The mixed tissue sample library was prepared using the GS Titanium Rapid Library Preparation Kit (Roche, Basel, Switzerland) following the manufacturer's protocol and was then analysed using a 454 Titanium platform, with a single end in one plate.

Data analysis.
For both RNA-Seq data and de novo assembled data, sequences were checked for the quality and were trimmed using the FASTX toolkit. De novo assembly was conducted using the Newbler programme 57 . Assembled contigs were BLAST queried and annotated to the A. thaliana database (27,416 protein data points in TAIR 10) using Blastx. The top hit gene with an e-value < = 1e−10, was treated as an ortholog. RNA-Seq reads were then mapped to the references and counted using Bowtie 58 . RNA-Seq reads were mapped to the references and the FPKM (fragments per kilobase of transcript per million fragments mapped) values were calculated. Because the chemicals used in Illumina sequencing differed between the first and second runs, the RNA-Seq data for plants (A and B) were analysed separately. FPKM values were normalized by the quantile-spline method using the normalize qspline function in the affy package of R 59 . Log 2 -transformed values of the normalized FPKM plus 2 -5 were used as expression values (Supplementary Data S1). Principal component analysis (PCA) was performed using the prcomp function in R to compare the overall gene expression patterns of the rhizome, shoot and root apexes and leaves. Genes whose expressions were higher than 0 in at least one sample were used in the PCA. We performed k-mean clustering using the k-mean function in R. The number of clusters was set to 16. Genes whose expressions were higher than 5 in at least one sample were used in the k-mean clustering. GO enrichment analyses were conducted for each cluster.
In order to characterize the transcriptomes of the four tissues of C. leucantha, especially that of the rhizome tip, the gene expression patterns were compared with previously reported transcriptome data from other species. First, comparisons were made with A. thaliana microarray data for seven selected tissues collected as part of the AtGenExpress project: the hypocotyl 7 days after germination (d.a.g.) (ATGE_2), the root at 7 d.a.g. (ATGE_3), the root at 17 d.a.g. (ATGE_9), the vegetative shoot apex at 7 d.a.g. (ATGE_6), the rosette leaf at 17 d.a.g. (ATGE_12), mature pollen 6 weeks after germination (ATGE_57) and the flower at stage 10/11 (ATGE_31) 60 . The Spearman's rank correlation between each of C. leucantha transcriptome and each A. thaliana transcriptome was calculated. Putative orthologous genes in A. thaliana were used in the calculation.
Anatomical observation and subcellular structure of rhizome tips. In order to anatomically characterize the rhizome tips of C. leucantha, bright-field microscopy and transmission electron microscopy were used. The rhizome tip samples were collected in June when the rhizome elongated actively from the transplants. The samples were fixed with 2% paraformaldehyde (PFA) and 2% glutaraldehyde (GA) in 0.05 M cacodylate buffer (pH 7.4) at 4 °C overnight. After fixation, the samples were rinsed three times with 0.05 M cacodylate buffer for 30 min each, followed by post fixation with 2% osmium tetroxide (OsO4) in 0.05 M cacodylate buffer at 4 °C for 3 h. The rhizome tips were observed using a transmission electron microscope (JEM-1200 EX; JEOL Ltd., Tokyo, Japan) at an acceleration voltage of 80 kV. Digital images were taken using a CCD camera (VELETA, Olympus Soft Imaging Solutions GmbH, Münster, Germany). To examine the presence or absence of ER bodies in leaves and rhizomes of C. leucantha, tissues were collected from plants transplanted in the CER garden (described above). The ER body is reported to be present in the roots and cotyledons and absent in the leaves and stems of A. thaliana 25 . A fluorescence producing gene, SP-GFP-HDEL, which is expected to localize GFP in the endoplasmic reticulum (ER) and ER-derived organelles 61,62 was transiently introduced into the collected tissues by particle bombardment, as described previously 63 . The tissues were then incubated for 2 days in a 22 °C chamber under continuous light with 100% humidity. GFP-expressing cells were inspected under a confocal laser-scanning microscope (LSM780 META; Carl Zeiss, Oberkochen, Germany).