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.
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 (Fig. 1a) and were subjected to RNA-Seq analysis. Hereafter, we refer to these two samples as plant A and plant B, respectively. Samples from the four tissues of plants A and B are referred to as rhizomes A and B, shoots A and B, leaves A and B, and roots A and B, respectively. These tissue samples were analysed using an Illumina HiSeq 2000 (Illumina, San Diego, CA, USA), using 2 × 101-bp paired-end sequencing, with one lane for two plants.
As a result of RNA-Seq, 10.9 M, 19.4 M, 17.7 M and 12.4 M reads in plant A and 21.8 M, 16.6 M, 15.2 M and 20.1 M reads in plant B were obtained for the rhizome, shoot, root and leaves, respectively. RNA-Seq data were mapped to the reference genes using Bowtie (version 0.12.8). For all samples, 65–80% of total reads were 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 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) were annotated to Pathogenesis-related protein 4 (PR-4, AT3G04720), Light-dependent Short Hypocotyls 10 (LSH10, AT2G42610), Peroxidase 37 (AtPrx37, AT4G08770), GDSL-like Lipase 23 (GLL23, AT1G54010), Peroxidase 39 (AtPrx39, AT4G11290), Methionine Sulfoxide Reductase B9 (MSRB9, AT4G21850) and β-glucosidase 23 (PYK10 /BGLU23, AT3G09260). PR-4 is involved in JA (Jasmonic Acid)-dependent defenses and is upregulated by treatment with necrotrophic fungi27 and rhizobacteria28. A homolog of LSH10 in potato (Solanum tuberosum L.) has been reported to show very high transcript levels in stolons (rhizomes) and young tubers29. 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 leaves30. 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 uptake31. 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 homeostasis32. MSRB9 has been reported to be prevalently expressed in A. thaliana roots and is involved in tolerance to the accumulation of ROS33. 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 alkanes34,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 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 seedlings36. In A. thaliana, ESP was found to be consistently present in the epidermal cells of all aerial parts37 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 Arabidopsis38. 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 Log2 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 expansion in the roots of A. thaliana39,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 light41. ARR3 has been reported to be constitutively expressed in the vascular tissue of both shoots and roots and is induced by cytokinin in root tissues42. 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 complex43,44, which is a huge protein aggregate that facilitates the enzymatic activity of PYK1045. 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 bodies48, 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 PYK1024, 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, but not chloroplasts (Fig. 4d), suggesting that the rhizome of C. leucantha serves as the energy storage organ49,50,51. 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.
Reproducibility of gene expression patterns
To evaluate the reproducibility of the results of RNA-Seq, we conducted real-time qPCR experiments using the three independent sets of four tissues, i.e. rhizomes, shoots, roots and leaves, from three plants. We examined ten selected homologous genes, i.e. LSH10 (isotig20431), Prx37 (isotig06471), GLL23 (isotig17255), PYK10 (isotig03131), BGLU18 (isotig00469), PBP1 (isotig17838), NAI2 (isotig07117), TSA1 (isotig02612), AGL9 (isotig04090) and PHOT1 (isotig12495). Gene expressions were reproduced in all genes at least in the rank order across four tissues except for those of PBP1 and PHOT1 (Fig. S1).
Expressions of additional eight genes., i.e. AGT (isotig04516), CA1 (isotig18828), CER (isotig00800), ER (isotig12468), EXT19 (isotig00056), MYB15 (isotig19642), PSBX (isotig27053) and PUB23 (isotig08161), whose patterns were extremely different between leaf (aboveground part) and root (belowground one) were also checked by real-time qPCR. The expression patterns were mostly similar to those of RNA-Seq (Fig. S2).