Skip to main content
Download PDF

Transcriptome analysis reveals plasticity in gene regulation due to environmental cues in Primula sikkimensis, a high altitude plant species

BMC Genomics volume 20, Article number: 989 (2019) Cite this article

Abstract

Background

Studying plasticity in gene expression in natural systems is crucial, for predicting and managing the effects of climate change on plant species. To understand the contribution of gene expression level variations to abiotic stress compensation in a Himalaya plant (Primula sikkimensis), we carried out a transplant experiment within (Ambient), and beyond (Below Ambient and Above Ambient) the altitudinal range limit of species. We sequenced nine transcriptomes (three each from each altitudinal range condition) using Illumina sequencing technology. We compared the fitness variation of transplants among three transplant conditions.

Results

A large number of significantly differentially expressed genes (DEGs) between below ambient versus ambient (109) and above ambient versus ambient (85) were identified. Transcripts involved in plant growth and development were mostly up-regulated in below ambient conditions. Transcripts involved in signalling, defence, and membrane transport were mostly up-regulated in above ambient condition. Pathway analysis revealed that most of the genes involved in metabolic processes, secondary metabolism, and flavonoid biosynthesis were differentially expressed in below ambient conditions, whereas most of the genes involved in photosynthesis and plant hormone signalling were differentially expressed in above ambient conditions. In addition, we observed higher reproductive fitness in transplant individuals at below ambient condition compared to above ambient conditions; contrary to what we expect from the cold adaptive P. sikkimensis plants.

Conclusions

We reveal P. sikkimensis’s capacity for rapid adaptation to climate change through transcriptome variation, which may facilitate the phenotypic plasticity observed in morphological and life history traits. The genes and pathways identified provide a genetic resource for understanding the temperature stress (both the hot and cold stress) tolerance mechanism of P. sikkimensis in their natural environment.

Background

Understanding constraints on species’ range limits have long been a primary goal of ecologists [1]. Climate has been recognized as a factor controlling species’ range limit [2]. When the climate changes gradually, ecosystems and species can evolve together. However, given the current rate at which the climate is changing [3], concerns are rising about the capacity of species to adapt. Sessile organisms such as plants have to be considerably more adaptable to stressful environments and must acquire greater tolerance to multiple stresses than animals. It is well known that environment induced phenotypic plasticity plays an important role in adaptation [4, 5], and plant phenotypic responses to altered environmental stresses are mainly regulated through gene expression [6, 7]. Thus, understanding plasticity in gene expression in natural systems is crucial, for predicting and managing the effects of climate change on plant species.

Variation in gene expression patterns plays a key role in the evolution of phenotypes [8] that allow an organism to acclimatize to stress [9, 10]. For example, thermal stress is considered a major constraint to plant reproduction. Almost all organisms respond to thermal stress by synthesizing heat-shock proteins (HSPs) [11,12,13]. However, different species respond differently to similar stress conditions; cold stress induces over expression of the C-repeat binding factor (CBF) genes in Arabidopsis thaliana [14] and induces over expression (10-fold upregulation) of OsCYP19–4 gene in Oryza sativa [15]. Plants may respond differently to multiple stress conditions [16], and the molecular mechanisms associated with multiple stresses might differ from those related to single stress [17, 18]. While many studies provide insight into plant responses to single stresses under controlled conditions [19,20,21], responses to changing conditions in the natural environment remains less understood.

Variation in gene expression under different conditions can be identified through genome-wide transcriptome analysis [22] using RNA sequencing (RNA_seq) [6, 23]. Application of RNA-seq to non-model species allows the use of their transcriptomes to understand their responses to changes in the environment [24, 25]. Many studies clearly demonstrated/ suggested that adaptive plasticity can processed through transcriptome variation [26,27,28,29], and much work is needed in these regards.

Altitudinal gradients provide a wide temperature range over a very short distance [30] and are therefore ideal to study potentially adaptive phenotypic variation in plants in the wild. Temperature differences along this fine-scale altitudinal gradients across ‘space’ can be used to infer the potential temporal responses of a population to climate change [31]. Many studies on altitudinal gradient to date have focused on species morphological and physiological differences, or the genetic basis of high altitude adaptations, and few studies have examined the contribution of gene expression level variation along altitudinal gradients [32, 26, 28]. Primula sikkimensis (genus Primula L.) is high altitude specialist plant, and one of the most dominant and widespread species, distributed along the altitudinal gradient of Sikkim Himalaya (27 °C 62’N, 88 °C 63’E) from 3355 m a.s.l. to 4598 m a.s.l. (field survey during 2012–2015, Lachen valley North-Sikkim). Populations sampled at different altitudes display phenotypic differences. Populations from higher altitudes (~ 4500 m a.s.l.) are smaller with delayed maturity and flowering compared to lower altitude populations (~ 3500 m a.s.l.), which are taller and flower earlier in the spring [33].

In this study we carried out a transplant experiments within and beyond the altitudinal range limit of P. sikkimensis. The gene expression profiles of transplant groups were obtained with transcriptome sequencing and we identified differentially expressed genes (DEGs) between within and beyond range transplant groups. The overall objective of this study was to facilitate a better understanding of how the gene expression level variation may have contributed to abiotic stress compensation in Primula sikkimensis.

Results

Illumina paired-end sequencing and de novo assembly of transcriptome

Illumina paired-end sequencing generated approximately 90 million raw reads (2 × 101 base pair). After pre-processing of raw reads, approximately 60 million reads (R1 = 2 × 94 base pair & R2 = 2 × 101 base pair) were left. In the absence of available reference genome for P. sikkimensis, we de novo assembled the transcriptome to be used as a reference for read mapping and gene expression profiling (hereafter referred to as the reference transcriptome assembly). We assembled the high-quality processed reads and the best-combined assembly resulted in 67,201 genes, 81,056 transcripts with a mean length of 785.87 bp and average open reading frame (ORF) length of 468.6 bp. The N50 of contigs was 1359 bp, a total size of 63.4 Mb, and a GC content of 38.99%. Similarly, results of separate assemblies in all the three transplant conditions were documented in Table 1. Only 3% (2647) of the transcripts have putative frameshifts which suggests good quality transcriptome data (Accession number: SRP150603). The raw reads generated from Illumina sequencing were deposited at National Centre for Biotechnology Information (NCBI), SRA with accession numberSRP150603.

Table 1 The results of separate transcriptome assemblies of P. sikkimensis in all three transplant conditions (ambient, below ambient and above ambient), and the reference assembly generated by combining the reads from all three conditions were documented in tabular form
Full size table

Functional annotation and identification of differentially expressed genes (DEGs)

Functional annotation of P. sikkimensis transcriptome assembly was carried out using TRAPID, in which Plaza database was used. Plaza is a collection of transcripts and genomes of plants. Our annotation resulted in 22,332 (27.6%) of transcripts annotated with GO categories and 26,313 (32.5%) of P. sikkimensis sequences annotated with known protein domains.

Using the RNA-seq data, we derived gene expression profiles in P. sikkimensis for all three transplant conditions. We then carried out two comparative transcriptome analyses between Ambient (A) the control, versus Below Ambient (BA), and Above Ambient (AA) transplant conditions. For comparison of differentially expressed genes we used 21,167 transcripts which mapped to the reference transcriptome of P. sikkimensis. To judge the significance of gene expression difference from the two pairwise comparisons we identified significantly differentially expressed genes of P. sikkimensis as those with log2 (fold change) ≥ 2 and log10 (p-value) < 0.05, as a threshold. A large fold change in expression does not always imply statistical significance, as those fold changes may have been observed in genes that received little sequencing or with many iso-forms [34], therefore we consider both fold change and p-value to identify the significant DEGs. We used volcano plots to show the significant DEGs which relate the observed differences in gene expression to the significance associated with those changes under Cuffdiff’s statistical model (Fig. 1). We found 109 significant DEGs from BA vs. A comparison, 81 up-regulated and 28 down-regulated (Fig. 2a).These genes include heat shock proteins HSP20, HSP70, Transcriptional factor B3, Methionine synthase, Zinc finger, dTDP-4-dehydrorhamnose reductase, DNA-binding, ATPase, and UDP-glucuronosyl (full list of genes, Additional file 8 Table S3a). From AA vs. A, we found 85 significant DEGs of which 61 were up-regulated and 24 were down-regulated (Fig. 2a). These genes include Heat shock protein DnaJ, bZIP transcription factor and Histone H5 (full list of genes, Additional file 8 Table S3b). Forty genes were common between the two pair-wise comparisons, whereas 69 and 45 genes were unique to BA vs. A and AA vs. A comparison respectively (Fig. 2b).

Fig. 1
figure 1

Volcano plots showing differentially expressed genes between (a) below ambient vs. ambient and (b) above ambient vs. ambient. The y-axis corresponds to the mean expression value of log10 (p-value), and the x-axis displays the log2 fold change value. The orange dots represent the significantly differentially expressed transcripts (p < 0.05); the black dots represent the transcripts whose expression levels did not reach statistical significance (p > 0.05)

Full size image
Fig. 2
figure 2

Differential gene expression profiles. a A number of up and down regulated genes in the pair-wise comparison between below ambient versus ambient and above ambient versus ambient transplant conditions. b Venn diagram presenting the number of unique and overlapping genes between two pair-wise comparisons

Full size image

Gene ontology (GO) and pathways mapping of DEGs

DEGs from the two pair-wise comparisons were mapped to GO database and GO terms were assigned. The DEGs had a GO ID and were categorized into small functional groups in three main categories (cellular component, molecular function, and biological process) of GO classification. Based on sequence homology, 42 and 36 functional groups were categorized in BA vs. A, and AA vs. A comparisons, respectively. Among these groups, “cell” and “cell part” were dominant within the “cellular component” category; “binding” and “catalytic” were dominant in the “molecular function” category; and “cellular process” and “metabolic process” were dominant in the “biological process” category (Additional file 4 Figure S4b).

The biological function associated with significant DEGs were further analyzed in terms of enriched Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways [35]. The DEGs had a KO ID and were categorized into small pathways. A total of 34 pathways were predicted for BA vs. A comparison and among them, “metabolic pathway”, “biosynthesis of secondary metabolites” and “flavonoid biosynthesis” were the most highly represented categories (Additional file 9 Table S4a). Similarly, 23 pathways were predicted for AA vs. A comparison and among them, “metabolic pathway”, “biosynthesis of secondary metabolites”, “plant hormones signal transduction”, and “photosynthesis” were the most highly represented categories (Additional file 9 Table S4b). The top 15 KEGG pathways of DEGs in these two pairwise comparisons are shown in Fig. 3.

Fig. 3
figure 3

Scatter plot of KEGG pathway enrichment analysis of differentially expressed genes in (a) below ambient versus ambient and (b) above ambient versus ambient transplant conditions. The number of DEGs in the pathway is indicated by the circle area, and the circle color represents the range of the corrected p-value (q-value) from 0~1. We display the top 15 pathway terms enriched by KEGG database

Full size image

Validation of RNA-Seq data by real-time quantitative RT-PCR

To confirm the RNA-Seq data, the transcript level of randomly selected 10 genes was examined by Real-Time quantitative RT-PCR (Fig. 4). All the genes exhibited the same pattern of expression as per FPKM (fragments per kilobase of exon per million fragments mapped) values for A, BA, and AA conditions except for “c15913_g1” annotated as ferredoxin-type protein, which was not detected in AA (Fig. 4). Taken together, all the selected genes (Table 2) showed same patterns that were consistent with the RNA-seq data, validating our experimental results.

Fig. 4
figure 4

Real-Time PCR analysis of selected genes in AA, A, and BA samples (a-j). Here the data repersented are realtive quantification (RQ) values of gene expression

Full size image
Table 2 List of primers used for Real-Time quantitative RT-PCR
Full size table

Differences in fitness-related traits of transplants across three transplant sites

Survival (rhizome sprouting) of transplants at the Ambient (A) the control site and Below Ambient (BA) transplant sites were > 85%, whereas the survival rate decreased to < 50% at Above Ambient (AA) site (Fig. 5a). We observed a significant decrease (Fig. 5b; ANOVA: F(2, 109) = 47.77, p < 0.001) in the height of P. sikkimensis outside of their range limit at BA and AA sites compared to A site. Post hoc comparisons using the Tukey HSDtest [36] indicates that the mean scores for the plant height at three transplant conditions was significantly different (BA: M = 22.41, SD = 10.96; A: M = 29.84, SD = 7.33; AA: M = 9.36, SD = 5.96). Similarly, flower number, representing the initial stage of reproductive fitness, also showed a significant decrease (Fig. 5c; ANOVA: F (2, 58) = 40.7, p < 0.001) outside the species range limit. Post hoc comparisons using the Tukey HSDtest [36] indicates that the mean scores for the flower number decrease significantly at BA and AA condition compared to A condition (BA: M = 6.08, SD = 2.92; A: M = 17.10, SD = 6.39; AA: M = 6.47, SD = 3.12). However, reproductive fitness represented by average seed production by transplants, was approximately seven seeds per individual at A and BA site, whereas the seed production dropped to four seeds per individual at AA site (Fig. 5d; ANOVA: F (2, 26) = 3.39, p = 0.05). Post hoc comparisons using the Tukey HSDtest [36] indicates that the mean scores for the seed production decreases significantly at AA (BA: M = 7.25, SD = 2.49; A: M = 7.50, SD = 3.00; AA: M = 4.66, SD = 2.12). Although seed production per individual was higher at A and BA site, the number of individuals producing seeds was less at BA site relative to A site. At A site 12 individuals produced seeds whereas at BA site only 8 individuals produced seeds. Similarly, at AA site, 9 individuals produced seeds. Taken together, we observed an overall decrease in fitness component of P. sikkimensis outside their present range limit (Fig. 4a-d), relative to range centre.

Fig. 5
figure 5

a Survival of transplanted rhizomes of P. sikkimensis at below ambient, ambient, and above ambient transplant sites. b plant height, c flower number and, d seed number: box plots showed differences among transplants at below ambient, ambient and above ambient transplant sites. Each box-and-whisker plot represents the observed measures for each population, with the centre bar indicating the median value. Bars with different letters are significantly different (Turkey post hoc tests, p < 0.05) and the numbers (n) above each bar of panel represents the sample size

Full size image

Discussion

Our gene expression analysis demonstrated that plastic gene expression variations have contributed to adaptation in high altitude Himalayan plant species (Primula sikkimensis) to different stresses in its natural environment. We identified a large number of genes with plastic expression differences between Ambient versus Below Ambient and Above Ambient conditions. The genes and pathways identified are good candidates for targeted studies of the role of variation in gene expression of a high altitude species to both the hot and cold temperature stress in its natural environment.

Are mechanisms of stress response conserved?

The below ambient and above ambient transplant sites are located beyond the altitudinal range limit of P. sikkimensis, with a temperature differences of approximately 2–3 °C (hotter) and approximately 1–6 °C (colder). Therefore, we compared the significant DEGs of P. sikkimensis from the BA vs. A comparison with heat stress genes of Arabidopsis thaliana using Gene Expression Omnibus (GEO), at National Center for Biotechnology Information (NCBI). Similarly, the genes from the AA vs. A comparisons were compared to the cold temperature stress genes of A. thaliana. Out of 109 significant DEGs of BA vs. A, 83 genes (76%) showed similarity with A. thaliana heat stress genes and out of the 85 genes from the AA vs. A comparison 56 genes (65.9%) were similar to A. thaliana cold stress genes (Thermal stress (hot): BA vs. A = 76% and (cold): AA vs. A = 65.9%). This supports previous work which suggests that the transcriptomic response to temperature stress might be highly conserved across plant species [37]. The plants at BA site with a higher temperature condition differentially up-regulated more genes than plants at AA site with a cold temperature condition; possibly indicating that expression of an elevated number of genes is necessary for the maintenance of P. sikkimensis individuals under heat stress conditions. This suggests that the high-temperature conditions, rather than the cold temperature conditions cause greater differences in the gene expression pattern of P. sikkimensis in our study.

How are below and above ambient different?

Plants are susceptible to adverse environmental conditions. Abiotic stresses such as extreme temperatures, drought, and high UV are some of the typical environmental stressors that can damage physiological functions, and reduce growth and yield of plants [38,39,40]. In plant communities, environmental stress can be a major source of plant mortality because plants are unable to escape from environmental stress through migration. Constant increases in ambient temperature are considered to be one of the most detrimental environmental stresses affecting plant growth and development [41]. Heat stress is not unique to plants and is also found in other organisms [42]. Heat stress at the molecular level causes an alteration in expression of genes involved in direct protection from high-temperature stress. These include genes responsible for the expression of osmoprotectants, detoxifying enzymes, transporters and regulatory proteins [13]. In our study, cytochrome P450, Pyridoxal phosphate-dependent decarboxylase, ubiquitin, transcriptional factor B3, HSPs, glycoside hydrolase family 16, NAD-dependent epimerase/dehydratase, haem peroxidize are some significant DEGs up-regulated in high-temperature conditions at BA transplant site. Similarly the cytochrome P450, Pyridoxal phosphate, ubiquitin, and glycoside hydrolase family are some of the genes which have been extensively studied in other plants in response to heat stress [43]. On the other hand Heat shock proteins (Hsp20, Hsp70), calcium-dependent protein kinase, glutamine aminotransferase are some significant DEGs down-regulated in high-temperature conditions at BA site (Fig. 1a). These results revealed that most of the genes involved in plant growth and development were up-regulated under BA conditions in P. sikkimensis whereas genes involved in signalling and stress-induced proteins (HSPs) were down-regulated. HSPs are proteins found in plant and animal cells in responsive to heat stress [44, 45]. HSPs generally functions as molecular chaperones, and are divided into HSP20, 40, 60, 70, 90, 100 and small HSP (sHSPs) [46]. HSPs have been shown to increase levels of gene expression when plants are exposed to elevated temperature [47]. However, our result revealed that HSP20 and HSP70 were down regulated by heat stress at BA site. As HSPs have been shown to be expressed more under heat stress over short time periods [48, 49] it seemed that in our study HSP20 and HSP70 genes might had responded for short time period after transplanting plants under heat stress at BA site but decreased with time.

Cold stress also adversely affects plant growth, development, and reproduction. Cold acclimation in plants involves reprogramming of gene expression [50]. Gene expression is induced by cold stress [51, 52] in a number of genes. These genes are thought to be involved in stress tolerance. In case of Arabidopsis, the protein kinases and transcription factors are some of the genes that are up-regulated in response to low temperatures [53]. In our study, Serine/threonine-protein kinase, phosphoinositide-binding, bifunctional inhibitor/plant lipid transfer protein/seed storage, transcription factor GRAS, DNA-binding WRKY are up-regulated in cold temperature conditions at AA site (Fig. 1b). These results revealed that most of the genes involved in signalling, defence and membrane transport/permeability were up-regulated under AA conditions in P. sikkimensis. It is evident from the pathway analysis that various genes involved in metabolic processes, secondary metabolism, and flavonoid biosynthesis were differentially expressed in BA, whereas genes involved in photosynthesis and plant hormones signalling were differentially expressed in AA site (Fig. 3).

Fitness variation within and beyond the range limit of P. sikkimensis

We observed a decline in rhizome sprouting of transplant individuals at the AA site (< 50% germination rate), whereas at A and BA sites rhizome sprouting was greater than 85%. The decreased in rhizome sprouting of P. sikkimensis at AA transplant site may be attributed to cold temperature, as the temperature is an important environmental variable regulating the sprouting of plant rhizomes [54]. Similar to our study, other studies also found the cold temperature to be an important factor in reducing and/or delaying the rhizome sprouting [54, 55]. For example, temperature range of 25–35 °C was optimum for the sprouting of turmeric (Curcuma longa L.) rhizomes, and sprouting did not occur below 10 °C [56]. We observed a decrease in plant height and flower number outside species range limit at BA and AA sites. However, the seed production (as a measure of reproductive fitness) was seven seeds per individual at both BA and A site but decreased to four seeds per individual at AA site. The higher reproductive fitness of the BA transplants compared to AA transplants was contrary to what we expect from the cold adaptive P. sikkimensis plants. According to the past distributional record of P. sikkimensis [57], the species was previously present at much lower elevation than its present lower elevational limits. Therefore, it is possible that the lower elevational limit of P. sikkimensis is not determined purely by abiotic factors, but that biotic factors are also be playing a role. These findings suggested that P. sikkimensis will favour hot temperature conditions for its germination and reproduction rather than cold conditions, if climate change proceeds in the high altitude of Sikkim Himalaya. However, expression of an elevated number of genes was necessary for the maintenance of P. sikkimensis individuals at the hotter temperature conditions as observed in case of BA transplant condition.

Conclusions

Plasticity in gene expression

Our study documented many differences in fitness related traits and gene expression associated with thermal stresses which suggest that P. sikkimensis undergoes a great deal of plasticity in its transcriptomic profiles. Transcriptomic plasticity of this species may facilitate the phenotypic plasticity in morphological and fitness related traits. Comparing the transcriptome profile of P. sikkimensis within and beyond the altitudinal range limit of species provided an opportunity to test for the plastic transcriptomic responses of species to stressful environmental condition specifically the thermal stresses. More importantly, transcriptome studies in naturally varying environments show that observed transcriptomic patterns may differ from those seen in controlled experimental conditions. Naturally varying environments may provide a better indication of responses of high elevation plants to ongoing climate change. However, despite the plastic responses of P. sikkimensis’ transcriptome, the transplant experiment resulted in reduced growth fitness and a decrease in the number of seed producing individuals under temperature stress conditions at BA and AA sites. This decrease indicates the vulnerability of species to future climate change. Future studies combining transcriptomic and genomic data may help in determining the evolutionary significance of transcriptomic variation responses to environmental stress and provide insights into plastic and evolutionary responses to climate change.

Methods

Transplant experiment0020

In order to conduct a direct test for adaptive significance of a phenotypic change [58, 59], we transplanted individuals of Primula sikkimensis in different environments. Our study was conducted along the elevation gradient of Lachen valley, Sikkim Himalaya (27 °C 62′N, 88 °C 63′E), India in three experimental sites (Additional file 1 Figure S1). The sites were classified as: ambient (A: 3951 (meters above sea level) m a.s.l.), below ambient (BA: 3256 m a.s.l.) and above ambient (AA: 4687 m a.s.l.) after validating elevation range of the species in field (Additional file 6 Table S1). Ambient site was within the altitudinal range limit of focal P. sikkimensis species and serves as a control for the experiment. The ambient site is located at an altitude of 3951 m a.s.l. locally knows as Thangu village with an average day temperature of 15 °C and average night temperature of 10 °C during the peak flowering seasons from June to August. Ambient site is also the source population of the P. sikkimensis rhizomes selected for the transplant experiment. Below ambient transplant site was approximately 100 m below the lower most elevational range limit of P. sikkimensis (3355 m a.s.l.), and the above ambient site was approximately 100 m above the upper elevational range limit of species (4598 m a.s.l.). There is a change in mean annual temperature of 1.5–3 °C at both the below ambient and above ambient transplant sites with respect to its lower most and uppermost elevational range limits, representing stressful environmental conditions (Additional file 2 Figure S2). According to the Intergovernmental Panel on Climate Change (IPCC) report the observed change in temperature from 1901 to 2012 in the Himalayas was 2–3 °C, therefore we designed our experiment to mimic this ongoing environmental change [60, 61]. A total of 300 rhizomes of P. sikkimensis were collected from the ambient site for the transplant experiment and were grown individually in a single pot of diameter 22.5 cm and height 30 cm. We kept 100 pots at the ambient site as a control for the experiment, and translocated 200 pots to the below and above ambient sites (100 at each site). We used soil from the ambient site to keep the variation in microfloral properties of the soil constant across transplant sites. The experiment was conducted from March 2013 to September 2014, and the fitness observations of transplants were made starting from March to September of 2014 (rhizome sprouting to reproductive phase). iButtons (hygrochron temperature data logger), were placed at each transplant site to measure the temperature at 2-h intervals from March to September 2014. Our temperature data showed a typical decreasing trend of temperature with increasing elevation as expected [62].

Transcriptome analysis

Plant material, RNA extraction, cDNA library synthesis and Illumina sequencing

Fresh leaf samples of P. sikkimensis growing at A, BA, and AA transplant sites were used for transcriptome studies. We collected leaf tissue samples in triplicates (one leaf sample * 5 individuals) randomly from each of the three transplant sites (A, AA, BA) in the same day and the samples were frozen in liquid nitrogen on-site and stored in − 80 °C till RNA isolation.

Total RNA was isolated from leaf samples collected in triplicates from three experimental conditions, using modified RNA isolation protocol [63]. RNA integrity was measured on 1% formaldehyde agarose gel by monitoring distinct 28S and 18S rRNA bands. Purity and concentration of isolated RNA were assessed by monitoring A260/A280 using NanoDrop spectrophotometer 2000C and Bioanalyzer (Additional file 7 Table S2). RNA samples (three RNA samples * three experimental conditions) with RNA Integrity Number (RIN) greater than eight were used for library preparation and sequencing. Transcriptome sequencing was performed using Illumina HiSeq1000 sequencing technology at the Next Generation Genomic Facility at the Centre for Cellular and Molecular Platforms (C-CAMP), Bangalore. The cDNA library for transcriptome sequencing was prepared using TruSeq RNA sample preparation kit V2 from Illumina, as per the manufacturer’s recommendations. The cDNA library was then sequenced using Paired-End 100 base pair chemistry using TruSeq PE cluster kit V3-cBot-HS and TruSeq SBS kit V3-HS for sequencing on the Illumina HiSeq 1000 platform following the manufacturer’s recommended protocols.

Read processing and de novo transcriptome assembly

FastQC and FASTX-Toolkit was used for quality checking and pre-processing of raw reads [61]. Read quality was checked and visualized with FastQC [64] and reads with a Phred scaled quality score of less than Q20 was removed. Reads were sorted using FASTX - Toolkit (Hannon Lab) and trimming of low-quality reads was performed using FASTX – Trimmer.

De novo transcriptome assembly of P. sikkimensis from all three conditions was performed independently by combining the filtered reads of biological triplicates. The reference transcriptome assembly was generated by combining reads of biological triplicates from all three conditions i.e., A, BA and AA. Denovo assembly was performed for all K-mers from 19 to 61 at an interval of two using Trinity [65]. Contigs shorter than 100 base pair (bp) were eliminated.

Transcript differential abundance calculation

Transcript abundance quantification was performed in terms of fragments per kilobase of exon per million fragments mapped (FPKM). A reference transcriptome (combining all three conditions) of non-redundant combined assembled transcriptome sequences at 90% sequence similarity by CD-HIT-EST [66] was taken. The reads from all the samples were mapped back to the reference transcriptome by using TopHat2 [67] at default parameters. Technical duplicates were merged using SAMTOOLS. Cufflinks was used to generate a GTF file for each gene model from the combined transcriptomic sequences [68]. The FPKM values for each transcript in all the samples were determined. The differential abundance of transcripts among different samples/sites was calculated by the cuffdiff 2 tool [69]. The complete workflow is provided in Additional file 3 Figure S3. Gene Ontology (GO) is an international standardized gene functional classification system which describes properties of genes and their products in any organism. GO is standardized gene functional classification system and it has three ontologies: cellular component, molecular function, and biological process. Functional annotation in terms of GO [70] and gene family of the transcripts was done by using TRAPID [71], an online server. TRAPID gives the option to search against the available databases viz., PLAZA 2.5 and OrthoMCLDB version 5. After getting GO annotation, the functional enrichment of the transcripts was also performed and abundant transcripts in each condition were plotted by using WEGO tool [72](Additional file 4 Figure S4).

Identification of differentially expressed genes (DEGs) and functional annotation

Analyses of DEGs includes the screening of genes that were differentially expressed among two pair-wise comparisons (BA vs. A and AA vs. A), and GO functional enrichment and KEGG (Kyoto Encyclopaedia of Genes and Genomes) pathway enrichment analysis for these DEGs. We used a value of log2 (fold change) ≥ 2 and mean expression value of log10 (p-value) < 0.05, as the threshold to judge the significance of gene expression difference [34]. We used Blast2GO [73] to get GO annotation for significant DEGs of two pair-wise comparisons. After getting GO annotation for every DEGs, we used WEGO [72] to do GO functional classification. KEGG is a pathway-related database and pathway enrichment analysis identifies significantly enriched pathways in DEGs [74], and KOBAS [75] was used to test the statistical significance of the enrichment of DEGs in KEGG pathways [35].

Real-time PCR analysis

Total RNA was isolated from A, BA, and AA samples as described above for transcriptome analysis. RNA extracts were treated with DNase I, amplification grade (Invitrogen, USA) to remove DNA contamination. Complementary DNAs (cDNAs) were synthesized using SuperScript III cDNA synthesis kit (Invitrogen, USA) as per the protocol.

The differentially expressed genes were selected randomly for qRT-PCR from A, BA, AA conditions based on their FPKM values (Additional file 5 Figure S5). The primers for all the genes were designed using the Primer3Plus software [76] as listed in Table 2. All the PCR reactions were performed in triplicatesin 10 μl reaction mixture containing diluted cDNA samples as template, 2× SYBR® Green Master Mix (Applied Biosystems, USA), and 200 nM each of forward and reverse gene-specific primers (Table 2). The reactions were performed in StepOnePlus™ Real-Time PCR System (Applied Biosystems, USA) using the following program: initial denaturation at 94 °C for 10 min, followed by 40 cycles of amplification (94 °C for 30s, 60 °C for 30s, and 72 °C for 30 s) and final melt curve analysis was performed. Transcript levels of all the genes were normalized with an internal control reference 26S rRNA gene [77]. The relative expression ratio of each gene was calculated using the comparative Ct value method as described previously [78]. Here, the transcript levels represented are relative quantitation (RQ) values of gene expression. Expression is shown after normalization to 26S rRNA gene. Values were calculated using the ΔΔCT method, and the error bars represented as RQMIN and RQMAX.

Survival, growth and reproductive fitness of transplant individuals among three transplant sites

We recorded the number of rhizomes sprouted from the transplanted rhizomes at each site. Height (cm) of transplants was measured as a representative of growth fitness. Flower number and seed number per transplanted individual in each transplant site were quantified as a measure of reproductive fitness. The differences in growth and reproductive fitness of transplant individuals between the three transplant sites (A, BA and AA) were assessed using ANOVA. Here we consider all the transplant individuals that manage to reach growth and/or reproductive maturity including the individuals which we used for RNA-seq analysis. When significant differences were observed, the ANOVA was followed by the Tukey posthoc tests [36] for pair-wise comparisons after Bonferroni correction.

Availability of data and materials

The data supporting the results of this article are available in the National Centre for Biotechnology Information (NCBI) [Accession number: SRP150603 (SRX4219916, SRX4219915, SRX4219914, SRX4219913, SRX4219912, SRX4219911, SRX4219910, SRX4219909, SRX4219908)]. All the supporting data are available in the electronic supplementary material.

Abbreviations

A:

Ambient

AA:

Above ambient

BA:

Below ambient

DEG:

Differentially expressed gene

GO:

Gene Ontology

KEGG:

The Kyoto Encyclopedia of Genes and Genomes

m a.s.l.:

Meters above sea level

bp:

Base pair

GEO:

Gene Expression Omnibus

NCBI:

National Centre for Biotechnology Information

FPKM:

Fragments per kilobase of exon per million fragments mapped

RIN:

RNA Integrity Number

IPCC:

Intergovernmental Panel on Climate Change

RQ:

Realtive quantification

References

  1. Brown JH, Stevens GC, Kaufman DM. The geographic range: size,shape, boundaries, and internal structure. Annu Rev Ecol Syst. 1996;27:597–623.

    Article  Google Scholar 

  2. Gaston KJ. Geographic range limits: achieving synthesis. Proc R Soc B Biol Sci. 2009;276:1395–406.

    Article  Google Scholar 

  3. Gaffney O, Steffen W. The Anthropocene equation. Anthr Rev. 2017;4:53–61.

    Google Scholar 

  4. Bradshaw AD. Evolutionary significance of phenotypic plasticity in plants. Adv Genet. 1965;13:115–55.

    Article  Google Scholar 

  5. Forsmsn A. Rethinking phenotypic plasticity and its consequences for individuals, populatins and species. Heredity. 2015;115:276–84.

    Article  Google Scholar 

  6. López-Maury L, Marguerat S, Bähler J. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nat Rev Genet. 2008;9:583–93.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  7. Hodgins-Davis A, Townsend JP. Evolving gene expression: from G to E to G ?? E. Trends Ecol Evol. 2009;24:649–58.

    PubMed  PubMed Central  Article  Google Scholar 

  8. Stern DL, Orgogozo V. The loci of evolution: how predictable is genetic evolution? Evolution. 2008;62:2155–77.

    PubMed  PubMed Central  Article  Google Scholar 

  9. Jaglo-Ottosen KR. Arabidopsis CBF1 overexpression induces COR genes and enhances freezing tolerance. Science. 1998;280:104–6.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  10. Xu C, et al. De novo and comparative transcriptome analysis of cultivated and wild spinach. Sci Rep. 2015;5:17706.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. Howarth BYCJ, Ougham HJ. Tansley Review No . 51 Gene expression under temperature stress. New Phytol. 1993;125:1–26.

    CAS  Article  Google Scholar 

  12. Swindell WR, Huebner M, Weber AP. Transcriptional profiling of Arabidopsis heat shock proteins and transcription factors reveals extensive overlap between heat and non-heat stress response pathways. BMC Genomics. 2007;8:125.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  13. Hasanuzzaman M, Nahar K, Alam MM, Roychowdhury R, Fujita M. Physiological, biochemical, and molecular mechanisms of heat stress tolerance in plants. Int J Mol Sci. 2013;14:9643–84.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  14. Cook D, Fowler S, Fiehn O, Thomashow MF. From the cover: a prominent role for the CBF cold response pathway in configuring the low-temperature metabolome of Arabidopsis. Proc Natl Acad Sci. 2004;101:15243–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  15. Yoon DH, et al. Overexpression of OsCYP19-4 increases tolerance to cold stress and enhances grain yield in rice (Oryza sativa). J Exp Bot. 2016;67:69–82.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  16. Nicky JA, Peter EU. The interaction of plant biotic and abiotic stresses: from genes to the field. J Exp Bot. 2012;63:3523–44.

    Article  CAS  Google Scholar 

  17. Rasmussen S, Barah P, Suarez-Rodriguez MC, Bressendorff S, Mundy J, et al. Transcriptome responses to combinations of stresses in Arabidopsis. Plant Physiol. 2013;161:1783–94.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. Sewelam N, Oshima Y, Mitsuda N, Ohme-Takagi M. A step towards understanding plant responses to multiple environmental stresses: a genome-wide study. Plant Cell Environ. 2014;37:2024–35.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  19. Li Y, Wang Y, Tang Y, Kakani VG, Mahalingam R. Transcriptome analysis of heat stress response in switchgrass (Panicum virgatum L.). BMC Plant Biol. 2013;13:153.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  20. Calzadilla PI, Maiale SJ, Ruiz OA, Escaray FJ. Transcriptome response mediated by cold stress in Lotus japonicus. Front Plant Sci. 2016;7:374.

    PubMed  PubMed Central  Article  Google Scholar 

  21. Chen S, Li H. Heat stress regulates the expression of genes at transcriptional and post-transcriptional levels, revealed by RNA-seq in Brachypodium distachyon. Front Plant Sci. 2017;7:1–13.

    Google Scholar 

  22. Nagalakshmi U, et al. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008;320:1344–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  23. Ekblom R, Galindo J. Applications of next generation sequencing in molecular ecology of non-model organisms. Heredity. 2011;107:1–15.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  24. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57063.

    Google Scholar 

  25. Alvarez M, Schrey AW, Richards CL. Ten years of transcriptomics in wild populations: what have we learned about their ecology and evolution? Mol Ecol. 2015;24:710–25.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  26. Cheviron ZA, Whitehead A, Brumfield RT. Transcriptomic variation and plasticity in rufous-collared sparrows (Zonotrichia capensis) along an altitudinal gradient. Mol Ecol. 2008;17:4556–69.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  27. Guevara DR, Champigny MJ, Tattersall A, Dedrick J, Wong CE, et al. Transcriptomic and metabolomic analysis of Yukon Thellungiellaplants grown in cabinets and their natural habitat show phenotypic plasticity. BMC Plant Biol. 2012;12:175.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. Yang W, Qi Y, Lu B, Qiao L, Wu Y, et al. Gene expression variations in high-altitude adaptation: a case study of the Asiatic toad (Bufo gargarizans). BMC Genet. 2017;18:62.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  29. Cho SM, Lee H, Jo H, Lee H, Kang Y, et al. Comparative transcriptome analysis of field- and chamber-grown samples of Colobanthus quitensis (Kunth) Bartl, an Antarctic flowering plant. Sci Rep. 2018;8:11049.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. Körner C. The use of ‘altitude’ in ecological research. Trends Ecol Evol. 2007;22:569–74.

    PubMed  Article  PubMed Central  Google Scholar 

  31. Hollister RD, et al. Corrections for Elmendorf et al., Experiment, monitoring, and gradient methods used to infer climate change effects on plant communities yield consistent patterns. Proc Natl Acad Sci. 2015;112:4156.

    Article  CAS  Google Scholar 

  32. Hovenden MJ, Vander Schoor JK. Nature vs nurture in the leaf morphology of southern beech, Nothofagus cunninghamii (Nothofagaceae). New Phytologist. 2003;161:585–94.

    Article  Google Scholar 

  33. Gurung PD, Ratnam J, Ramakrishnan U. Facilitative interactions among co-flowering Primula species mediated by pollinator sharing. Plant Ecol. 2018;219:1159–68.

    Article  Google Scholar 

  34. Trapnell C, et al. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2012;31:46–53.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  35. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  36. Tukey JW. Comparing individual means in the analysis of variance. Biometrics. 1949;5:99–114.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  37. Zhang B, Pan X, Cannon CH, Cobb GP, Anderson TA. Conservation and divergence of plant microRNA genes. Plant J. 2006;46:243–59.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  38. Boyer JS. Plant productivity and environment. Science. 1982;218:443–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  39. Xiong L, Schumaker KS, Zhu JK. Cell signaling during cold, drought, and salt stress. Plant Cell. 2002;14:S165–83.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  40. Krasensky J. Jonak C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory.

  41. Lobell DB, Field CB. Global scale climate–crop yield relationships and the impacts of recent warming. Environ Res Lett. 2007;2:14002.

    Article  Google Scholar 

  42. Schlesinger MJ, Ashburner M, Tissieres A. Heat shock : from Bacteria to man: Cold Spring Harb. Lab; 1982.

  43. Zhang X, et al. Expression profile in Rice panicle: insights into heat response mechanism at reproductive stage. PLoS One. 2012;7:e49652.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. Ritossa F. A new puffing pattern induced by temperature shock and DNP in drosophila. Experientia. 1962;18:571–3.

    CAS  Article  Google Scholar 

  45. Jung K, Gho H, Nguyen MX. Genome-wide expression analysis of HSP70 family genes in rice and identification of a cytosolic HSP70 gene highly induced under heat stress. Funct Integr Genomics. 2013;13:391–402.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  46. Waters ER. The evolution, function , structure , and expression of the plant sHSPs. J Exp Bot. 2013;64:391–403.

    CAS  Article  Google Scholar 

  47. Aitken SN, et al. Hierarchical responses of plant-soil interactions to climate change: consequences for the global carbon cycle. Ecol Lett. 2015;8:127–33.

    Google Scholar 

  48. Silva-Correia J, Freitas S, Tavares RM, Lino-Neto T, Azevedo H. Phenotypic analysis of the Arabidopsis heat stress response during germination and early seedling development. Plant Methods. 2014;10:1–11.

    Article  Google Scholar 

  49. Song K, Yim WC, Lee B. Expression of Heat Shock Proteins by heat stress in Soybean. Plant Breed Biotech. 2017;5:344–53.

    Article  Google Scholar 

  50. Viswanathan C, Zhu J. Molecular genetic analysis of cold-regulated gene transcription. Philos Trans R Soc Lond B Biol Sci. 2002;357:877–86.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. Laroche Â, Frick M, Huel Â, Puchalski B, Gaudet DA. Cold induced expression of plant defensin and lipid transfer protein transcripts in winter wheat. Physiol Plant. 2003;117:195–205.

    Article  Google Scholar 

  52. Janska A, Marsik P, Zelenkova S, Ovesna J. Cold stress and acclimation – what is important for metabolic adjustment? Plant Biol. 2010;12:395–405.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  53. Knight MR, Knight H. Low-temperature percpetion leading to gene expression and cold tolerance in higher plants. New Phytol. 2012;195:737–51.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  54. Loddo D, Masin R, Otto S, Zanin G. Estimation of base temperature for Sorghum halepense rhizome sprouting. Weed Res. 2012;52:42–9.

    Article  Google Scholar 

  55. Satorre EH, Rizzo FA, Arias SP. The effect of temperature on sprouting and early establishment of Cynodon dactylon. Weed Res. 1996;36:431–40.

    Article  Google Scholar 

  56. Yukio I, Md AH, Keiji M, Hikaru A, Takuji H. Effects of planting date on emergence, growth and yield of turmeric (Gurcuma longa L.) in Okinawa prefecture, southern Japan. Jpn J Trop Agr. 2004;48:10–6.

    Google Scholar 

  57. Polunin O, Stainton A. Flowers of the Himalaya: Oxford University press; 1984.

  58. Hereford JA. Quantitative survey of local adaptation and fitness trade-offs. Am Nat. 2009;173:579–88.

    PubMed  Article  PubMed Central  Google Scholar 

  59. Blanquart F, Kaltz O, Nuismer SL, Gandon S. A practical guide to measuring local adaptation. Ecol Lett. 2013;16:1195–205.

    PubMed  Article  PubMed Central  Google Scholar 

  60. IPCC (Intergovernmental Panel on Climate Change). In: Solomon S, et al., editors. Climate change 2007: the physical science basis. Contribution of Working Group I to the Fourth assessment report of the Intergovernmental Panel on Climate Change. Cambridge: Cambridge University Press; 2007. p. 235–336.

    Chapter  Google Scholar 

  61. IPCC. In: Stocker TF, et al., editors. Climate Change 2013: The Physical Science Basis Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. New York: Cambridge University Press; 2013. p. 1535.

    Google Scholar 

  62. Gale J. Availability of carbon dioxide for photosynthesis at high altitudes: theoretical considerations. Ecology. 1972;53:494–7.

    Article  Google Scholar 

  63. Muoki RC, Paul A, Kumari A, Singh K, Kumar S. An improved protocol for the isolation of RNA from roots of tea (Camellia sinensis (L.) O. Kuntze). Mol Biotechnol. 2012;52:82–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  64. Andrews S. Babraham Bioinformatics - FastQC A Quality Control tool for High Throughput Sequence Data. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

  65. Grabherr MG, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  66. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  67. Kim D, et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:36.

    Article  CAS  Google Scholar 

  68. Trapnell C, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  69. Trapnell C, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7:562–78.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  70. Ashburner M, et al. Gene ontologie: tool for the unification of biology. Nat Genet. 2000;25:25–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. Van Bel M, et al. TRAPID: an efficient online tool for the functional and comparative analysis of de novo RNA-Seq transcriptomes. Genome Biol. 2013;14:134.

    Article  Google Scholar 

  72. Ye J, et al. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34:293–7.

    Article  Google Scholar 

  73. Conesa A, et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  74. Ogata H, et al. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 1999;27:29–34.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  75. Wu J, Mao X, Cai T, Luo J, Wei L. KOBAS server: a web-based platform for automated annotation and pathway identification. Nucleic Acids Res. 2006;34:720–4.

    Article  CAS  Google Scholar 

  76. Untergasser A, Nijveen H, Rao X, Bisseling T, Geurts R, Leunissen JA. Primer3Plus, an enhanced web interface to Primer3. Nucleic Acids Res. 2007;35(Web Server issue):W71–4. https://doi.org/10.1093/nar/gkm306.

    Article  PubMed  PubMed Central  Google Scholar 

  77. Singh K, et al. 26S rRNA-based internal control gene primer pair for reverse transcription-polymerase chain reaction-based quantitative expression studies in diverse plant species. Anal Biochem. 2004;335:330–3.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  78. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods. 2001;25:402–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgments

Thanks to Forests, Environment & Wildlife Management Department, Government of Sikkim; Department of Home, Government of Sikkim; Army and Police department of Sikkim for the research and road permits in the Kanchenjunga Biosphere Reserve, North-Sikkim, India. We also thank the field assistants and people of Lachen valley, North-Sikkim for assisting in the field and logistic support. Thanks to Bhusan Gurung and Ophilia Ibapalei Lyngdoh for technical assistance. Sincere thanks to Megan Aylward for her diligent proofreading of this manuscript.

Funding

Financial support for this study including fellowship and field work was supported by Department of Biotechnology (DBT), Government of India under the project “Technological Innovations and Ecological Research for the Sustainable use of Bioresources in the Sikkim Himalaya”. The funding body played no role in the design of the study and collection, analysis and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

  1. National Center for Biological Sciences (NCBS), Tata Institute of Fundamental Research, GKVK Campus, Bellary Road, Bengaluru, Karnataka, 560065, India

    Priya Darshini Gurung, Atul Kumar Upadhyay, Ramanathan Sowdhamini & Uma Ramakrishnan

  2. Manipal University, Manipal, India

    Priya Darshini Gurung

  3. Present Address: Thapar Institute of Engineering & Technology, Department of Biotechnology, Patiala, Punjab, 147004, India

    Atul Kumar Upadhyay

  4. Institute of Bioresource & Sustainable Development, A National Institute under Department of Biotechnology, Ministry of Science & Technology, Government of India, Gangtok, Sikkim, 737102, India

    Pardeep Kumar Bhardwaj

  5. Present address: Institute of Bioresources and Sustainable Development, Meghalaya, 6th Mile, Upper Shillong, Meghalaya, 793009, India

    Pardeep Kumar Bhardwaj

Authors
  1. Priya Darshini Gurung
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Atul Kumar Upadhyay
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Pardeep Kumar Bhardwaj
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Ramanathan Sowdhamini
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Uma Ramakrishnan
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

UR and PDG design the study. PDG performed the field experiment, sampling and primer designing for RT-PCR. PKB performed the RNA isolation, and RT-PCR. PDG, AKU, and RS analysed the data. PDG initiated the manuscript. All authors commented on the manuscript and gave approval for accountable of the final version of the manuscript.

Corresponding author

Correspondence to Priya Darshini Gurung.

Ethics declarations

Ethics approval and consent to participate

All fieldwork and plant sample collection were conducted legally. This study does not involve any species at risk of extinction. Field experiment and plant sample collection were approved by Forests, Environment & Wildlife Management Department, Government of Sikkim; and Department of Home, Government of Sikkim.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1.

Transplant sites.

Additional file 2: Figure S2.

Temperature across transplant sites.

Additional file 3: Figure S3.

Workflow of differential expression of genes.

Additional file 4: Figure S4.

WEGO.

Additional file 5: Figure S5.

Heatmap plot for 25 random DEGs.

Additional file 6: Table S1.

Details of three transplant sites.

Additional file 7: Table S2.

Yield and quality of RNA samples.

Additional file 8: Table S3.

Details list of significant DEGs in Below ambient versus Ambient and Above ambient versus Ambient transplant condition.

Additional file 9: Table S4.

Pathway details.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Gurung, P.D., Upadhyay, A.K., Bhardwaj, P.K. et al. Transcriptome analysis reveals plasticity in gene regulation due to environmental cues in Primula sikkimensis, a high altitude plant species. BMC Genomics 20, 989 (2019). https://doi.org/10.1186/s12864-019-6354-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-019-6354-1

Keywords

  • Gene expression
  • Transplant experiment
  • Transcriptomics
  • Climate change
  • Range limits

BMC Genomics

ISSN: 1471-2164

Contact us