Skip to main content
Download PDF

A comparison of heat-stress transcriptome changes between wild-type Arabidopsis pollen and a heat-sensitive mutant harboring a knockout of cyclic nucleotide-gated cation channel 16 (cngc16)

BMC Genomics volume 19, Article number: 549 (2018) Cite this article

Abstract

Background

In flowering plants, the male gametophyte (pollen) is one of the most vulnerable cells to temperature stress. In Arabidopsis thaliana, a pollen-specific Cyclic Nucleotide-Gated cation Channel 16 (cngc16), is required for plant reproduction under temperature-stress conditions. Plants harboring a cncg16 knockout are nearly sterile under conditions of hot days and cold nights. To understand the underlying cause, RNA-Seq was used to compare the pollen transcriptomes of wild type (WT) and cngc16 under normal and heat stress (HS) conditions.

Results

Here we show that a heat-stress response (HSR) in WT pollen resulted in 2102 statistically significant transcriptome changes (≥ 2-fold changes with adjusted p-value ≤0.01), representing approximately 15% of 14,226 quantified transcripts. Of these changes, 89 corresponded to transcription factors, with 27 showing a preferential expression in pollen over seedling tissues. In contrast to WT, cngc16 pollen showed 1.9-fold more HS-dependent changes (3936 total, with 2776 differences between WT and cngc16). In a quantitative direct comparison between WT and cngc16 transcriptomes, the number of statistically significant differences increased from 21 pre-existing differences under normal conditions to 192 differences under HS. Of the 20 HS-dependent changes in WT that were most different in cngc16, half corresponded to genes encoding proteins predicted to impact cell wall features or membrane dynamics.

Conclusions

Results here define an extensive HS-dependent reprogramming of approximately 15% of the WT pollen transcriptome, and identify at least 27 transcription factor changes that could provide unique contributions to a pollen HSR. The number of statistically significant transcriptome differences between WT and cngc16 increased by more than 9-fold under HS, with most of the largest magnitude changes having the potential to specifically impact cell walls or membrane dynamics, and thereby potentiate cngc16 pollen to be hypersensitive to HS. However, HS-hypersensitivity could also be caused by the extensive number of differences throughout the transcriptome having a cumulative effect on multiple cellular pathways required for tip growth and fertilization. Regardless, results here support a model in which a functional HS-dependent reprogramming of the pollen transcriptome requires a specific calcium-permeable Cyclic Nucleotide-Gated cation Channel, CNGC16.

Background

Temperature stress is a major contributor to crop loss around the world, with pollen infertility being one of the most important underlying causes [1,2,3,4,5,6]. Fertilization during plant reproduction is highly sensitive to hot and cold temperatures, with even a single hot day or cold night carrying the potential to disrupt reproductive success. Given the uncertainties of climate change, understanding this vulnerability is significant because most of the world’s food supply is derived from seed crops that depend on pollen fertilization.

Plant cells are proposed to sense heat stress (HS) through several mechanisms [7,8,9,10], including 1) accumulation of ROS (reactive oxygen species), 2) temperature-responsive calcium (Ca2+) channels, 3) an endoplasmic reticulum (ER)-unfolded protein response (UPR), 4) increased membrane fluidity, and 5) increased transcriptional activities related to a temperature-sensitive binding affinity of histones to specific regions of chromatin, leading to access of transcriptional regulators that can alter transcription.

Several transcriptome studies have been conducted to gain insights into heat stress responses (HSRs) in flowering plants during reproductive development, including both microarray and RNA-Seq experiments (reviewed in [11, 12]). However, only a few studies examined isolated mature pollen grains, and no pollen HS-studies have been reported for Arabidopsis. Nevertheless, studies in rice [13,14,15] and tomato [16,17,18,19] support an expectation that pollen and other reproductive tissues respond to HS with global changes in gene expression, including changes in the abundance of transcripts encoding small heat shock proteins (HSPs) [14,15,16, 18], heat shock transcription factors (HSFs, e.g. HsfA2) [14, 16, 18], enzymes involved in ROS mitigation (such as ascorbate peroxidase 1 (APX1) and catalase 2 (CAT2)) [14,15,16], proteins associated with the ER-unfolded protein response (such as HSPs, e.g. Hsp90 and CDC48, a homohexameric AAA-ATPase chaperone) [16, 18], and membrane trafficking (such as vesicle-associated membrane proteins AtVAMP725 and AtVAMP726) [16]. In addition, RNA-Seq has allowed the discovery of HS-dependent changes in small non-coding RNAs (sncRNAs) and alternative splicing events in tomato pollen [17, 18].

A previous study from Tunc-Ozdemir et al. [20] identified Arabidopsis Cyclic Nucleotide-Gated cation Channel 16 (CNGC16) as a necessary component for pollen fertility under multiple stress conditions, including HS. A loss-of-function mutation of cngc16 (i.e., knockout) resulted in a 10-fold reduction in pollen fitness and seed production under a HS condition. To obtain mechanistic insights into this hypersensitivity, we performed an RNA-Seq experiment to compare the pollen transcriptomes from a cngc16 mutant and WT Col-0, with and without a temperature stress (hot/cold stress regime described in [20]). Results here define an extensive HS-dependent reprogramming of approximately 15% of the WT pollen transcriptome, and identify at least 27 transcription factor associated changes that could provide unique contributions to a pollen HSR. In contrast to WT, pollen from a cngc16 mutant showed 1.9-fold more HS-dependent transcriptome changes. These results support a model in which HS-tolerance in pollen involves an extensive reprogramming of the transcriptome through a signaling pathway that is dependent on the function of a specific Cyclic Nucleotide-Gated cation Channel, CNGC16.

Results

To compare cngc16 and WT pollen for differences in their response to a temperature stress condition, transcriptomes were analyzed from mature pollen grains harvested at midday from plants grown under control (normal) conditions or a HS regime. For the stress condition, plants were grown under a diurnal cycle of hot and cold temperatures (Additional file 1) [20] for a period of 1 week, and pollen were harvested at the end of the HS period that peaked at 40 °C. Thus, many of the harvested pollen grains matured under multiple cycles of hot and cold stress conditions. As a result of pollen developing under this continuous stress cycle, their transcriptomes are predicted to reflect an adaptation to and “memory” of these conditions, as well as responses during the last 1 h HS period. While the mature pollen grains from Arabidopsis are desiccated, it is not clear how responsive they are to transcriptome changes during this “dormant state”, and to what degree the HS-changes observed here reflect the final stages of the stress regime, or the history of the stress.

A total of 12 pollen samples were processed with three independent biological replicates for each genotype and growth condition. Transcriptome sequencing was done with a single Illumina flow cell. Expression read counts (Additional file 2) were normalized (Additional file 3) to facilitate a comparison of all 12 samples. More than 21 million reads were obtained for each sample (Additional file 4a), with a principal component analysis (PCA) of the filtered and normalized expression data showing that 87% of the observed variance of the samples can be explained by differences in the stress states. (Principal Component 1, Additional file 4b). The expressed transcript fragments (reads) were aligned to the Araport 11 reference genome [21] for all non-obsolete gene models. While expression of at least three read counts were detected for 24,860 genes (Additional file 2), these results were filtered to identify 14,226 genes showing expression levels deemed high enough for quantification of transcripts in the following 8 gene-type categories annotated by Araport [21]: protein_coding, ara11_novel genes, long_noncoding_RNA, antisense_long_noncoding_RNA, miRNA, other_RNA, small_nuclear_RNA, small_nucleolar_RNA.

Two approaches were used to evaluate purity of pollen samples used for transcriptome comparisons. The first was to use light microscopy to visually inspect pollen samples for non-pollen debris that might have passed through a 70 μm nylon mesh screen during the harvesting protocol. This visual inspection indicated that filtered samples were free of any large fragments of vegetative tissues or non-pollen debris (data not shown).

In the second method (Additional file 5), transcriptomes were evaluated for the relative abundance of selected transcripts, and then compared to previously published pollen transcriptome studies, including both RNA-Seq and microarray data sets [22, 23]. These comparisons were done using two normalization strategies, each providing equivalent results. Data sets for comparison were normalized using either 1) total transcriptome read counts, or 2) a set of four housekeeping genes that appeared to show highly consistent expression levels across all 12 samples in this current transcriptome study (Additional file 5). The relative transcript abundances were evaluated for i) a control group of genes comprising three members of the CNGCs (CNGC7, 8, and 18), which were previously studied in the Harper lab and known to be essential to pollen fertility and show moderate to low levels of expression [24, 25], and ii) five known photosynthetic marker genes, three corresponding to light harvesting complex genes, and two corresponding to chlorophyll a/b-binding protein genes [26]. Although pollen grains are not considered to be photosynthetic, they contain plastids, with photosynthetic genes showing low basal levels of expression. To use these reference genes to evaluate whether the current study had a level of pollen purity similar to two previous pollen studies [22, 23], we confirmed that the average abundance ratio for three “pollen-specific” CNGC reference genes ranged from 1.0 to 2.1 (values of samples here divided by normalized values in [22] or [23]), which suggests that the pollen purity of the current study was equivalent (“1.0”) or even 2-fold more pure (“2.1”). As further support, the transcript abundance for the five different photosynthetic marker genes were all very low, and ranged from being 3-fold less to 1.7-fold higher in the current study compared to [22, 23], respectively. These comparisons indicate that the purity of the pollen in the current study is similar to that of other studies.

To test whether samples from a cngc16 knockout showed an expected deficiency in full-length CNGC16 mRNA, individual RNA-Seq reads corresponding to CNGC16 transcripts were aligned with a reference genome sequence (Additional file 6) using a Web-based tool, Integrative Genome Browser (IGB, bioviz.org). This alignment failed to detect any evidence of full-length transcripts in pollen harboring a cngc16–2 transfer DNA (T-DNA) insertion. Together, the absence of full-length transcripts along with a transgene rescue experiment from Tunc-Ozdemir et al. [20] supports the contention that cngc16–2 is a true loss-of-function knockout (i.e., null allele).

WT pollen HSR showed 2102 changes

To identify genes involved in a WT pollen HSR, we compared normalized transcript abundance profiles for pollen harvested from plants grown with or without a temperature stress regime described in [20]. This comparison showed 2102 statistically significant changes (≥ 2-fold changes with adjusted p-value ≤0.01; Fig. 1a; Additional file 3). To confirm the reliability of the RNA-Seq data, four genes were chosen and subjected to a quantitative PCR (Q-PCR) analysis (Additional file 7a, with RNA-Seq results also shown in Table 1). The Q-PCR analysis confirmed relative transcript abundance differences with an overall Spearman Correlation Coefficient computed as 0.72. As expected, correlations were highest for individual examples in which RNA-Seq results showed ≥2-fold differences. When comparing the WT-heat and cngc16-heat responses individually, the Spearman Correlation Coefficients were 1.0 and 0.67, respectively. The correlation dropped to 0.02 for the cngc16 control (normal). We believe this comparison was negatively impacted by the Myb29 transcripts being undetectable in cngc16 pollen under control (normal) condition, but still capable of being amplified into a detectable signal using Q-PCR.

Fig. 1
figure 1

A comparison of WT and cngc16 transcriptomes reveals numerous differences in their HSRs. a. Table shows number of significant changes in transcript abundance out of 14,226 quantifiable transcripts in WT and cngc16 pollen. b. Venn diagram showing number of significant HS-dependent changes that are shared or specific to WT or cngc16 pollen. c. A Venn diagram showing number of significant HS-dependent changes in transcription factors (TFs) that are shared or specific to WT or cngc16 pollen. All numbers are extracted from Additional file 3. Separate lists of transcript changes between unstressed WT and cngc16 mutant (21), specific and common HSRs (471, 2305, and 1631), and direct HSR (192) are provided in Additional file 9. Specific lists of TF changes are provided in Additional file 11

Full size image
Table 1 Transcriptome profiles for 33 transcription factor genes with potential relevance to a pollen HSR
Full size table

Given the ability of RNA-Seq strategies to detect and accurately quantify transcripts with very low abundance, an important question to be addressed on a gene-by-gene basis is whether a transcript with low-abundance has biological importance. It is certainly possible that significant biological impacts during normal development or HS can be associated with very small changes in transcript abundance, or even for transcripts that are below the detection limit of the current study. To give a perspective here on transcripts considered to be “low abundance”, the depth of sequencing in this study allowed for a minimum relative read count after normalization of 0.43 for the WT control (normal) samples. For these samples, the median abundance of all transcripts quantified was 31, or approximately 72-fold higher than the minimum threshold read count. One of the transcripts verified here by Q-PCR was chosen because it corresponded to a relatively low abundance of 1.5 in WT control (normal) conditions and showed a significant HS-dependent increase in abundance in the cngc16 mutant (i.e., transcription factor gene MYB29, Additional file 7a, with RNA-Seq results also shown in Table 1). This example provides confidence that transcripts showing changes near the lower limits of detection in this RNA-Seq data set can still be quantified with statistical confidence. Nevertheless, reliable detection does not address whether a low abundance transcript has an important biological function, or is simply present at near-background levels from “leaky” expression. As evidence here to emphasize the potential biological importance of relatively low abundance transcripts in pollen, we identified two examples of genes, with genetic evidence for a biological function in pollen, and abundance levels near or below our threshold for quantification. For example, at a low expression average of 7.2 in WT control samples, there is a gene encoding a putative acetyl-transferase for which we have genetic evidence for a biological role in pollen HS tolerance, based on a reduced pollen transmission efficiency and gene rescues for two independent T-DNA gene disruptions (Harper unpublished). In addition, despite an abundance level that failed to even rise above our threshold for quantification (e.g., 0.43 in WT), there is genetic evidence for a critical pollen autonomous function for AT4G01220 (MGP4), which encodes a sugar nucleotide transferase [27]. Thus, in the context of speculating on gene functions, these two examples provide a precedent that genes with relatively low abundance levels of 7.2 or lower can still have significant biological impacts on pollen fertility.

To evaluate whether the large number of transcriptome changes observed for the WT HSR could be linked to changes with transcription factors, we identified 89 HS-dependent changes associated with transcription factors in WT pollen (i.e., 23 + 66 in Fig. 1c), corresponding to 14% of the 616 transcription factors that were deemed quantifiable in our pollen transcriptomes. Of those, there were 27 showing at least a two-fold higher expression in pollen compared to seedling tissues (Table 1, based on ratios calculated from [22]). These pollen biased transcription factors represent potential regulatory nodes of importance to understanding the pollen-specific features of a HSR.

Heat stress also resulted in an increased abundance of 79 non-protein coding transcripts in WT (Table 2 and Additional file 3) including six microRNAs: MIR156A, MIR156C, MIR160, MIR836A, MIR831A, and MIR780A. Candidate targets for these microRNAs were determined by psRNATarget (http://plantgrn.noble.org/psRNATarget/, [28]) (see Additional file 8) and included members of the Squamosa Promoter Binding Protein-Like (SPL) and Auxin Response Factor (ARF) transcription factor families, which are known to contribute to stress tolerance in vegetative tissues [29] and regulate floral development [30]. Table 2 shows four categories of non-coding RNAs, with examples illustrating the two highest x-fold changes (up and down) for each category. Given that our RNA sample preparations included a selection step for poly-adenylated transcripts, there was an expectation that many of the non-coding RNAs would be selected against and therefore under-represented in our RNA-Seq results. This was corroborated by the observation that from a list of 689 different tRNAs (transcribed by RNA polymerase III), there were only four different transcripts detected with abundance levels that met our stringency thresholds for quantification (Additional file 2). Nevertheless, polyadenylation does occur for some non-coding RNAs that are transcribed by RNA polymerase II [31]. For example, a transcript encoding a miRNA might accumulate as a polyadenylated mRNA before processing into a mature miRNA, and these maturation events could be sensitive to heat stress. Thus, it is possible that some of the non-coding mRNAs identified in Table 2 include transcripts that accumulate before a final processing event. Regardless, their consideration as potential HS-dependent changes are warranted, as the statistical criteria imposed here still provide confidence that changes in their relative transcript abundances were quantified with a reasonable degree of reliability (p ≤ 0.01).

Table 2 Examples of HS-dependent changes in non-coding RNAs
Full size table

The cngc16 pollen HSR includes 3936 transcriptome changes

To define the transcriptome changes associated with a HSR in cngc16 pollen, we compared normalized transcript abundance profiles from pollen harvested from mutant plants grown with or without a temperature stress regime in parallel with WT plants discussed above. In cngc16 pollen, there were 3936 HS-dependent changes (≥ 2-fold changes with adjusted p-value ≤0.01; Fig. 1a; Additional file 3), which represents a 1.9-fold increase in the number of changes compared to the 2102 changes that met the same stringency criteria in WT.

Under control conditions there were only 21 significant differences between WT and cngc16 transcriptomes

To determine if the cngc16 knockout mutation resulted in a transcriptome with significant pre-existing differences under control conditions (i.e., “pre-existing condition”), we compared WT and cngc16 transcript profiles for pollen harvested from plants grown under normal conditions. Only 21 transcripts were significantly different based on the standard threshold criteria used here (≥ 2-fold changes with adjusted p-value ≤0.01; Fig. 1a; Additional file 9a). Nevertheless, these 21 examples still included large log2-fold differences that ranged from 3.7 to − 4.2. While the number of pre-existing differences is relatively small, it remains possible that one or more of these 21 differences potentiates a different HSR in the cngc16 mutant.

Under HS-conditions, there were 2776 differences between the WT and cngc16 HS responses

In contrast to the small number of transcriptome differences (i.e., 21) between WT and cngc16 pollen under normal conditions, there were 2776 differences in a simple contrast comparison between the two lists of genes independently categorized as HS-dependent in WT or cngc16 (Fig. 1b; Additional file 9b and c). These included 471 transcript changes that were unique to WT (Additional file 9b), and 2305 that were unique to cngc16 (Additional file 9c). However, to identify the most statistically significant differences, a direct comparison was made between the abundance of transcripts in the WT and cngc16 HS-transcriptomes. Using the same standard criteria applied above for normal conditions (≥ 2-fold changes with adjusted p-value ≤0.01), the HS-transcriptomes showed 192 differences (Additional file 9e, or Additional file 3 column R), of which 10 were already present under normal conditions (i.e., 182 new HS-dependent differences). This represents an approximate 9-fold increase in the differences between WT and cngc16 transcriptomes in response to a HS. Using a less stringent cut-off criteria that allowed for any difference between WT and cngc16 with a p-value greater than 0.05, there was a much greater number of 1531 differences between the two HS-transcriptomes, or an approximate 13-fold increase compared to control (normal) conditions. Thus, while there are potentially as many differences as there are similarities between the WT and cngc16 HSRs, a smaller subgroup of 182 HS-dependent differences can be categorized with a greater degree of statistically confidence (see Additional file 3, column R or Additional file 9e).

To highlight some of the most dramatic differences, Table 3 shows the top 23 examples of HS-dependent changes in WT that displayed the largest magnitude differences between WT and cngc16 under HS. These differences were analyzed for potential functional associations using STRING [32], revealing a cluster of seven genes that were classified as under-responsive to HS in cngc16 (Table 3; Additional file 10). While the functional significance of this under-responsive cluster is not known, it includes 4 genes related to cell walls and membrane dynamics. Interestingly, cell walls and membrane dynamics are the most frequent functional category for all of the genes in Table 3 (half of the top 20 largest differences).

Table 3 HS-dependent changes in WT that were most different in cngc16 under HS
Full size table

The WT and cngc16 HSRs showed 121 differences corresponding to transcription factors

A focused comparison of HS-dependent transcript changes for transcription factors was done to assess whether differential expression of transcription factors might contribute to the greater number of overall transcript changes in cngc16. While 66 HS-dependent transcription factor changes were common to cngc16 and WT pollen, there were 121 (1.8-fold more) changes categorized as potential differences (Fig. 1c; Additional file 11). However, in a direct comparison between the HS-transcriptomes of WT and cngc16, there were only seven transcription factor differences that satisfied a more stringent statistically criteria of ≥2-fold difference with an adjusted p-value ≤0.01 (see Table 1). Only one of these seven showed a significant HS-dependent change in both WT and cngc16 pollen (AT2G20110). This transcription factor also showed preferential expression in pollen compared to seedlings (listed in middle section of Table 1). The other six transcripts either just showed HS-dependent changes in cngc16 pollen, or were simply significant differences that were independent of a stress condition (bottom of Table 1). Regardless, any small changes in the abundance of transcription factors between WT and cngc16 could easily cause the larger number of downstream transcriptome differences observed under a HSR.

HS-dependent transcriptome changes for Ca2+-signaling related genes

Because CNGC16 is thought to function in both creating and sensing Ca2+-signals, a transcriptome comparison was done to evaluate whether WT and cngc16 pollen showed significant HS-dependent changes for a subset of genes associated with Ca2+-signaling (Ca2+-signaling related genes). Pollen expression was detected for 230 Ca2+-signaling related genes (Additional file 12). Of those, 40 (or 17%) showed significant changes in a WT HSR (Fig. 2). In a direct comparison between the WT and cngc16 HS-transcriptomes, there were no examples (except CNGC16) of any differences that met our standard criteria of ≥2-fold change with an adjusted p-value ≤0.01. This suggests that in the context of feedback networks that regulate Ca2+-signaling related genes, the loss of CNGC16 did not have a significant impact.

Fig. 2
figure 2

HS-dependent changes for 25 Ca2+-signaling related genes. Graph shows selected examples of HS-dependent log2-fold abundance changes for transcript encoded by 25 genes related to Ca2+-signaling. The 25 examples were selected as the most statistically significant changes in WT pollen (open bar) based on adjusted p-values. Corresponding transcript abundances from cngc16 mutant (grey bar) and aerial parts of seedlings (black bar) are shown for comparison. HS-dependent changes in transcript abundance in seedlings (aerial parts) was based on [33]. Complete list of Ca2+-related genes with putative functions (total 230) is provided in Additional file 12. CHXs are cation/H+-exchangers. CBLDs are Ca2+-dependent lipid-binding (CaLB domain) family protein. SLP is Calcineurin-like metallo-phosphoesterase superfamily protein. CRTs are calreticulins. ACA is autoinhibited Ca2+-ATPase pump. CAMs are calmodulins. CML is calmodulin like. CNGC is cyclic nucleotide gated channel. CPKs are Ca2+-dependent protein kinases. CBN is Ca2+-binding endonuclease/exonuclease/phosphatase family. CBL is calcineurin B subunit-like protein. CIPK is calcineurin B-like (CBL)-interacting protein kinase. CAX is cation exchanger. * indicates absence of probeset in microarray experiment

Full size image

To evaluate whether Ca2+-signaling related genes showed similar or different HSRs between pollen and vegetative tissues, a comparison was made with publicly available data sets obtained from AtGenExpress (http://jsp.weigelworld.org/expviz/expviz.jsp; [33]). While caution is required in making comparisons between experiments done with different HS conditions and transcript profiling technologies, the HS-dependent changes in pollen and seedlings appeared very different. For the subset of 40 Ca2+-signaling related genes that showed a significant change in WT pollen, only two genes showed an analogous HSR in seedling tissues (Calcium-dependent lipid-binding (CaLB domain) family protein, AT1G53590 and Calcium-binding endonuclease/exonuclease/phosphatase family, AT5G54130) (Fig. 2). Of the remaining 38, half did not even show a transcript abundance change in the same up or down direction. Thus, there was only a 5% overlap between HS-dependent changes in Ca2+-related genes detected in pollen and seedlings.

Discussion

While thermotolerance in plants has been widely studied, relatively little is known about specific effects of temperature stress on pollen [7, 11, 12]. Nevertheless, fertilization during plant reproduction is highly sensitive to hot and cold temperatures, and pollen is often considered a weak-link in reproductive stress tolerance [1,2,3,4,5,6].

This study provides a reference data set that identifies at least 2102 transcriptome changes associated with a HSR in WT Arabidopsis pollen (≥ 2-fold changes with adjusted p-value ≤0.01; Fig. 1a; Additional file 3). In addition, a parallel comparison with pollen from a cngc16 knockout mutant provides evidence that the HS-sensitivity caused by a cngc16 mutation [20] could be due in part to a defect in reprogramming the pollen transcriptome during a HSR.

HSRs in pollen and vegetative tissues have important similarities and differences

A potentially important conceptual difference between stress-tolerance programs in pollen and vegetative cells is that vegetative cells can often respond to a HS by down-regulating metabolism into a “survival mode”, and thereby wait for better growth conditions, whereas pollen tubes are often under temporal constraints in which there is a limited window of time to find and fertilize ovules. Because pollen must grow as fast as possible to successfully compete and “win a race” to fertilize a limited number of ovules, for some pollen types, a stress-response is postulated to include important differences compared to other cell-types in which the simplest solution to a HS is to down-regulate metabolism and wait for more optimal growth conditions.

The 2102 HS-dependent transcriptome changes observed here for WT pollen represent a dramatic reprogramming of at least 15% of the transcriptome (at ≥2-fold with adjusted p-value ≤0.01). The full extent of this reprogramming is most likely higher, with an estimate of 24% using a lower threshold criteria that includes smaller magnitude changes (i.e., less than two-fold), and/or abundance changes with a slightly more permissive p-value (e.g., ≤ 0.05 instead of ≤0.01). This extensive reprogramming has been observed in other studies in both plant and non-plant systems [34,35,36,37,38].

Many of the expected HS-dependent changes were observed here in pollen, such as an increased abundance of transcripts corresponding to heat shock transcription factors (HSFs, e.g. HsfA2), small heat shock proteins (HSPs), BCL-2-associated athanogene 6 (BAG6), WRKY transcription factors, Multiprotein Bridging Factor 1c (Mbf1C), dehydrogenases, phospholipases, and hormone-responsive transcription factors involved in ethylene (AT5G47230) and auxin (AT3G23050) signaling [20, 39,40,41,42,43]. Among the common HSRs were also antioxidant enzymes such as ascorbate peroxidase 2, peroxidase, and catalase [44,45,46].

Despite similarities with other HSRs, the following four observations create an expectation of important differences between a HSR in pollen compared to other plant cells. First, of the 89 transcription factors that showed a HS-dependent change in WT pollen, 84 (or 94%) failed to either be detected or show a similar HS-response in aerial parts of Arabidopsis seedlings (Additional file 11), based on a comparison with publicly available microarray data in AtGenExpress (http://jsp.weigelworld.org/expviz/expviz.jsp; [33]).

Second, pollen failed to show a HS-dependent change for 96% of a group of 67 genes that were curated as multi-stress regulated genes based on a comparison of nine abiotic stresses, including HS in Arabidopsis root and shoots [47]. Of these 67 genes, transcripts from only 19 (less than half) were detected in pollen (Additional file 13). Of those 19, only three showed a significant HS-dependent change in pollen, and of those three, two actually showed an opposite change in abundance compared to a vegetative HSR. The only example that showed a parallel HS-dependent decrease was AT1G05400, which is annotated as encoding a hypothetical protein.

A third observation was the failure to see expected changes in abundance of transcripts encoding SPL transcription factors (Squamosa Promoter Binding Protein-Like). In vegetative tissues, there is evidence that a HS-dependent increase in MIR156 triggers a decrease in transcripts for multiple SPL transcription factors [29, 48]. While the pollen HSR here also showed a HS-dependent increase in transcripts harboring MIR156a and c (seven to four-fold, respectively), a corresponding decrease in target abundance was not observed for any of the expected SPL genes (Additional file 8). In a second similar example, pollen showed a HS-dependent increase in transcripts harboring MIRNA160, but without a significant decrease detected in predicted regulatory targets (Additional file 8). These examples suggest that MIR156 and MIR160A might function differently during a HSR in pollen compared to other tissues.

A fourth observation was the poor correspondence between pollen and seedlings for HS-dependent responses among a subset of 230 pollen-expressed genes that are implicated in sensing or creating Ca2+-signals (Additional file 12). Of 40 Ca2+-signaling related genes with transcripts showing a HS-dependent change in WT pollen, there were only two examples of a similar response in vegetative tissues (Fig. 2). Likewise, in a subset of 23 HS-dependent changes in WT that were most different in cngc16 (Table 3), there was only one that was a potential HS-gene in seedlings, AT5G25280, a gene with unknown functions. Thus, together these four transcriptome-based observations above provide strong evidence that a HSR in pollen has significant differences compared to other cell types in plants.

A HSR in cngc16 pollen was significantly different than WT

Results here provide evidence that the HSR in cngc16 pollen is significantly different than WT. First, in a simple tally of genes categorized as HS-dependent in WT and cngc16, there were 2776 differences, with cngc16 pollen showing 1.9-fold more differences compared to WT (Fig. 1a; Additional file 3). Based on a more rigorous statistical analysis, a subset of 192 genes were found to have a ≥ 2-fold difference (with adjusted p-value ≤0.01) in a direct comparison of the HS-transcriptomes from WT and cngc16 (Additional file 3, column R and Additional file 9e).

While it is not yet clear which transcriptome differences between WT and cngc16 might have biological importance, some examples from the subgroup of the 192 most statistically significant differences include the following. First, there were seven transcription factors, of which one failed to be detected in the cngc16 transcriptome under HS or control (normal) conditions (AT2G34440, AGAMOUS-like 29). Because transcription factors regulate the expression of other genes, any transcription factor differences between WT and cngc16 could potentially amplify differences throughout the transcriptome during a HSR.

Second, in evaluating the 23 HS-dependent changes in WT that were most different in cngc16 (Table 3 and Additional file 10), the two most frequent functional categories were cell wall and membrane dynamics. In the context of cell wall, FAR3 (Fatty Acid Reductase 3, AT4G33790) showed increases in transcript abundance in cngc16 pollen that were opposite to WT. While FAR3 is involved in cuticular wax production in leaves [49], and is essential for pollen fertility [50], its potential importance to a HSR in pollen has not been investigated.

Another noteworthy example in the cell wall category was a gene encoding Leucine-Rich Repeat/Extensin 9 (LRX9, AT1G49490,). Studies on loss-of-function knockouts for members of this extensin subfamily have recently provided evidence for a critical role of these proteins in pollen tube growth [51,52,53,54]. Thus, a HS disruption in the expression of these, and other cell wall related genes, could potentially affect the structural dynamics of the cell wall and the ability of pollen tubes to grow and fertilize ovules.

In the context of membrane dynamics, there were three examples of genes encoding lipid transfer proteins (LTPs). While the specific biochemical activities of these proteins and their importance to a HSR is not understood, there is evidence that overexpression of an LTP in potato can protect cell membrane integrity under multiple stress conditions, including heat [55]. In addition, there are many studies showing the importance of regulating lipid biogenesis pathways for cells to adapt to temperature stresses [56]. As a specific example for pollen, a double knockout of lipid fippases ala6 and ala7 in Arabidopsis results in a stress-dependent pollen sterility [57] under the same temperature stress regime used here.

Role of CNGC16 in HSR

There are two general models to explain the function of CNGC16 in mediating a normal HSR. The first is a direct role of CNGC16 as part of an ion signaling pathway that functions in sensing or responding to HS. For example, a HS-triggered increase in cyclic-nucleotide monophosphate (cNMP) concentrations could activate CNGC16 and generate a cytosolic Ca2+ signal. While CNGC16’s ion selectivity properties have not been corroborated, a Ca2+ conductance has been suggested based on analogy to other closely related homologs [20]. For example, there is electrophysiology evidence for Ca2+ conductance from a study on a closely related cngc6 knockout [58]. In this case, CNGC6-mediated Ca2+-transients are thought to be critical to establishing vegetative HS-tolerance. In addition, CNGCa and CNGCd in a moss Physcomitrela have been implicated in modulating Ca2+ signals during a HSR [59]. Furthermore, a trio of isoforms closely related to CNGC16 in Medicago truncatula (CNGC15a, b, and c) have been implicated in controlling an elicitor-induced Ca2+ oscillation associated with the nucleus [60]. These Medicago homologs appear to be localized to the nuclear envelope, raising the potential that a HS-dependent Ca2+ signal in pollen might also be associated with the ER or nuclear envelope instead of the plasma membrane. However, regardless of subcellular location, a potential CNGC16-mediated Ca2+ signal could alter numerous cellular enzymes and structures, including the activity of key transcription factors important to a HSR.

A second alternative model is that a loss of CNGC16 might result in a cell with a pre-existing condition in which mutant cells are “unhealthy” and therefore less-able to sense or respond to multiple stresses, including a HS. While the transcriptome analysis here identified at least 21 differences that might contribute to an unhealthy pre-existing condition, it is not yet clear if any of these differences are actually responsible for the HS hypersensitivity of cngc16 pollen.

Regardless of a direct or indirect mechanism to explain why a cngc16 mutation results in pollen with hypersensitivity to HS, the HS-transcriptome response in cngc16 was clearly different than WT (e.g., Table 3, and Additional file 3, column R). In addition to the large number of differences, the nature of those differences suggests several possible mechanistic explanations for hypersensitivity, as revealed by comparing transcript changes in WT and cngc16 in the context of Gene Ontology (GO-term) classifications (e.g., Biological Process using PANTHER [61]) (Fig. 3 and Additional file 14). First, the cngc16 mutant showed a 1.8 to 2.2-fold increase in the number of HS-dependent changes categorized as responses to oxidative stress or abiotic stimulus. This increase is consistent with a model in which cngc16 pollen might be less-capable than WT in mitigating a ROS increase, which often occurs under stress situations. Second, cngc16 pollen showed an 8-fold greater number of differences in “regulation of cell growth”, and a 1.8-fold greater number of genes related to pollen development. These examples are consistent with a model in which cngc16 pollen are less-able than WT during a HSR to stabilize transcript levels associated with critical developmental programs related to cell growth in general, and pollen development in particular. While the cngc16 HSR compared to WT showed differences in more than 80% of 2380 different GO categories (Additional file 14), Fig. 3 includes two examples of controls in which transcript profiles for WT and cngc16 were similar (e.g., ATP metabolic processes and exocytic processes). A separate GO analysis (Additional file 15) was also done for the subset of 192 genes identified here showing the greatest differences in a direct comparison of cngc16 and WT under HS (Additional file 9e, or Additional file 3 column R). Under Biological Process, more than two thirds of the over-represented genes belonged to just two categories, metabolic process (GO:0008152) and cellular process (GO:0009987). In a subset of the 20 HS-dependent changes in WT that were most different in cngc16 (Table 3 and Additional file 10), the two most common categories were cell wall and membrane dynamics.

Fig. 3
figure 3

A GO (Gene Ontology) classification reveals similarities and differences between HS-triggered transcriptome changes in WT and cngc16 pollen. The GO analysis for Biological Process using PANTHER [61] was conducted for 2102 and 3936 HS-dependent changes in WT and cngc16 mutant (Additional file 14). The number of annotated genes from the input list that were mapped into a particular category is shown above each bar. Brackets show X-fold difference between cngc16 and WT. The categories shown are exocytic process (GO:0140029), ATP metabolic process (GO:0046034), response to oxidative stress (GO:0006979), response to abiotic stimulus (GO:0009628), and pollen development (GO:0009555). * and ** indicate enrichment in the GO category with p-value ≤0.05 and ≤ 0.01, respectively. A more comprehensive comparison of 5053 different GO categories is provided in Additional file 14

Full size image

Thus, HS appears to differentially perturb a large number of the transcriptome homeostasis networks in cngc16 pollen, with potential functional implications for a wide range of cellular and metabolic processes, including response to oxidative stress, regulation of cell growth, membrane dynamics and cell wall integrity. In this context, the simplest model to explain the HS-hypersensitivity of cngc16 pollen is that mutant cells are simply more vulnerable to losing control over multiple cellular systems, of which a failure with one or some combination make it impossible to carry out a very complicated development program.

Conclusions

The HS-transcriptome comparisons here provide new insights into a temperature-stress response in WT pollen. While pollen exhibit many of the same HS-dependent changes characteristic of vegetative tissues, evidence here suggests many important differences. For example, of the 89 transcription factors that showed a HS-dependent transcript abundance change in WT pollen, 94% failed to show a similar HSR in aerial parts of Arabidopsis seedlings (Additional file 11). Importantly, 27 of these transcription factor genes are expressed primarily in pollen (Table 1), and therefore represent regulatory nodes of potential importance to understanding pollen-specific features of a HSR.

A comparison of pollen from WT and a heat-sensitive cngc16 mutant showed only a small number of differences under control (normal) conditions (21), with differences increasing by at least 9-fold under HS. Given that more than 15% of the pollen transcriptome showed changes in response to a HS, it will be difficult to determine which specific differences in cngc16 pollen represent a potential cause of hypersensitivity, or simply arose as symptoms of a mutant cell’s inability to cope with HS. Nevertheless, results here suggest that cngc16 pollen have a defect in enacting or maintaining an appropriate HS-transcriptome response. Together, these results support a model in which reprogramming the pollen transcriptome for HS-tolerance is dependent on the proper functioning of a specific Cyclic Nucleotide-Gated cation Channel, CNGC16.

Methods

Plant material and growth conditions

The cngc16–2 (SAIL_726_B04, seed stock no. 91) mutant carries a glufosinate (Basta) resistance marker in the T-DNA insertion [20]. Following sterilization, seeds from both Arabidopsis Columbia WT (Col-0) and the cngc16–2 mutant were sown on 0.5× Murashige and Skoog (MS) medium (pH = 5.7) containing 1% agar with or without Basta (10 μg ml− 1). After 48 h of stratification at 4 °C, seedlings were grown at room temperature under constant light for 24 h for 10 days. The resulting seedlings were then transferred to soil (Sunshine SMB-238; Hummert). Unless otherwise stated, seedlings were grown until maturity under non-stress (control/normal) conditions in a growth chamber (Percival Scientific, Inc., http://www.percival-scientific.com/) with a photoperiod of 16 h of light and 8 h of dark at 22 °C, 40% humidity, and light intensity of 125 μmol m− 2 s− 1.

Stress condition

The stress regime with hot days and cold nights was performed as described in [20] and diagramed here in Additional file 1. Briefly, flowering plants grown under non-stress (control/normal) condition were transferred to a hot and cold stress-regime chamber and grown for one week before harvesting pollen. Pollen were harvested at the end of a 1 h HS-period at 40 °C (noon) following a gradual temperature rise from − 1 °C starting at 6 AM.

Pollen collection

For both the WT and mutant, three biological replicates per condition were used. Both HS and control (normal) pollen samples were harvested at midday. To purify pollen, open flowers were cut, vortexed for 10 to 20 s in a 50 ml centrifuge tube containing water, and filtered through a 70 μm nylon mesh (Becton Dickinson and Company). The pollen passing through the filter was centrifuged into a pellet for 30 s at 14000 rpm. The supernatant was discarded, and pollen pellets were frozen in liquid nitrogen and stored in − 80 °C.

RNA isolation, library preparation, and sequencing

Total RNA was extracted from pollen samples using Plant Mini Kit (Qiagen, Invitrogen), including an optional cleaning step using RNase-free DNase to eliminate genomic DNA contamination. One microgram of total RNA from each sample was sent on dry ice to the UCLA Neuroscience Genomic Core (UNGC) (Los Angeles, CA, USA) for library preparation using the TruSeq RNA kit (Illumina, San Diego, U.S.A) and sequencing. The RNA samples were quantified using the RiboGreen assay (ThermoFisher Scientific), and the quality of the RNA was checked by the Agilent TapeStation 2200. For each library, one microgram of total RNA was poly-A selected using oligo-dT magnetic beads, fragmented, and cDNA synthesized with random primers. Double-stranded cDNA was phosphorylated, and A-tailed followed by adapter ligation, PCR amplification, and sequencing. The cDNA libraries were multiplexed and run on a single lane. Paired-end (PE) sequencing was performed on an Illumina HighSeq2500 with two separate runs, the first generating sequences of 2 × 50 bp length and the second generating 2 × 69 bp read pairs.

Both ends of paired sequences were trimmed using Trimmomatic, version 0.36 [62] to remove sequences containing Illumina sequencing adapters, low-confidence bases (phred Q < 5), and sequences with length < 35 nucleotides. Sequence quality was measured and visualized before and after trimming using FastQC, version 0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). To verify the reverse/forward strand orientation expected of read pairs produced by TruSeq Stranded mRNA library preparation, subsets of 10,000 read pairs were randomly selected from trimmed read files with seqtk, version 1.0-r82 (https://github.com/lh3/seqtk), and compared to Araport 11 cDNA sequences [21] with BLASTN search via the command-line ncbi-blast+ (v2.5.0) application [63]. The plus/minus strand orientations of the tabular results were used to verify library forward/reverse orientations to correctly configure alignment tools.

The filtered sequence pairs were aligned to the Araport 11 version [21] of the TAIR10 Arabidopsis thaliana reference genome sequence [64], which was indexed in combination with coordinate information for positions of all exons and splice sites of transcribed genes annotated in the Araport 11 reference set, version 1.10.4 (released 06/2016) [21], using the spliced-read alignment tool HISAT2, version 2.0.5 [65]. Alignments produced by HISAT2 were converted to the binary BAM format and sorted with samtools, version 1.3.1 [66]. From these alignments, the number of read pairs aligned to known Arabidopsis genes was calculated using the FeatureCounts tool within the subread package, version 1.5.1 [67]. Alignments were counted once per pair, summarized at gene loci features, and read pairs with reported alignments to multiple loci were excluded from count totals. With at least three read counts set as a minimum threshold-limit for detection, there were 24,860 transcripts identified with correspondence to gene models or genomic features, excluding any matches to obsolete loci (Additional file 2).

To identify a subset of detectable transcripts deemed appropriate for quantification, two filtering steps were used. First, transcripts were excluded if they were categorized as rRNAs, tRNAs, transposons, or pseudogenes. The remaining categories annotated in Araport were retained for further processing: protein coding, ara11_novel genes, long_noncoding_RNA, antisense_long_noncoding_RNA, miRNA, other_RNA, small_nuclear_RNA, and small_nucleolar_RNA. To eliminate transcripts with expression levels considered too low for reliable quantification, transcripts with less than 10 fragments (counts) observed in at least two biological replicates in any condition were excluded. This left 14,226 transcripts for downstream analyses (Additional file 3). Data were normalized using the standard median ratio method for RNA-Seq data [68]. Principal component analysis (PCA) was performed on the normalized and filtered zero-centered counts per million data using singular value decomposition to validate clear separation between the different conditions (Additional file 4).

Differential gene expression

Differential gene expression between the four conditions was examined using DESeq2 [68]. Four comparisons including WT_heat vs. WT_control (normal), cngc16_heat vs. cngc16_control (normal), cngc16_control (normal) vs WT_control (normal), and cngc16_heat vs WT_heat were considered using simple contrasts. A multiple testing correction was performed within each of the four comparisons to adjust for the false discovery rate [69]. Genes with ≥2-fold changes and adjusted p-value ≤0.01 were considered to meet the standard significance threshold for this study.

Validation of RNA-Seq data by quantitative PCR (Q-PCR)

To validate RNA-Seq data, the transcript levels of four genes (Additional file 7) were examined by quantitative PCR (Q-PCR). Three of these four genes were chosen because they showed significant HS-dependent changes for the WT and or cngc16 samples in the RNA-seq analysis. The fourth gene was expressed at very low levels with a level of variation that made it a non-significant change. The same RNA samples used for RNA-Seq were used for verification of selected transcript changes using Q-PCR. First strand cDNA was synthesized using one microgram of total RNA via iScript cDNA Synthesis Kit (Catalogue#170–8891; Bio-Rad laboratory). A fraction (0.14 μg) of the cDNA was used as template in a 20 μL Q-PCR reaction mixture using SsoFast Evagreen Supermix (Catalog#172–5201; Bio-Rad laboratory) following the manufacturer’s instructions. Primer sequences used are shown in Additional file 7. Two reference genes cyclin p2;1 (CYCP2;1: AT3G21870) and UBIQUITIN 7 (UBQ7: AT2G35635) were chosen based on their minimal variation in the RNA-Seq analysis among all 12 experimental samples tested. The four other genes were selected to test examples of different patterns of changes observed in a comparison of WT and cngc16 transcriptomes. These included two transcription factors NAC105 (AT5G66300) and MYB29 (AT5G07690), as well as lipid transfer protein 6 (LTP6; AT3G08770) and delta vacuolar processing enzyme (DELTA-VPE; AT3G20210). Transcript abundance was quantified by Q-PCR (CFX96; Bio-Rad laboratory) with a separate normalization to the two different reference genes. The Q-PCR conditions were as follows: 30 s at 95 °C for enzyme activation, 39 cycles of 95 °C for 10 s for denaturing, and 25 s at 60 °C for annealing/extension.

A fold-change was calculated for each gene (normalized separately to each of the two reference genes UBQ7 and CYCP2) in relation to the expression of the WT control (normal) using the 2-ΔΔCT method [70] (Additional file 7). Based on all conditions and comparisons, the Spearman Correlation Coefficient between the Q-PCR and RNA-Seq expression values was computed as 0.72.

Gene ontology (GO)

Differentially expressed transcripts showing ≥2-fold changes and adjusted p-values ≤0.01 were analyzed using PANTHER [61]. Specifically, a statistical overrepresentation test (release 2017–04-13) was performed with the GO Biological Process Complete Annotation Data Set and a Bonferroni correction for multiple testing. The PANTHER Version 12.0 (release 2017–07-10) and Gene Ontology (GO) database (release 2017–08-14) were used. The test was based on the Arabidopsis genome of 27,060 annotated genes (released 07/2016).

Association network analysis

The top 23 examples of HS-dependent changes in WT with the largest magnitude differences between WT and cngc16 under HS were analyzed for potential associations or interactions using STRING version 10.5 [32] and MapMan version 3.6.0.[71].

Abbreviations

ACA:

Autoinhibited Ca2+-ATPase

APX:

Ascorbate peroxidase

ARF:

Auxin Response Factor

Avg:

Average

BAG6:

BCL-2-associated athanogene 6

Ca2+ :

Calcium

CaLB:

Calcium-dependent lipid-binding

CAM:

Calmodulin

CAT:

Catalase

CAX:

Cation exchanger

CBL:

Calcineurin B subunit-like protein

CBLD:

Ca2+-dependent lipid-binding

CBN:

Ca2+-binding endonuclease/exonuclease/phosphatase family

CHX:

Cation/H+-exchanger

CIPK:

Calcineurin B-like (CBL)-interacting protein kinase

CML:

Calmodulin like

CNGC:

Cyclic nucleotide-gated cation channel

cNMP:

Cyclic-nucleotide monophosphate

Col-0 WT:

Columbia-0 wild type

CPK:

Ca2+-dependent protein kinases

CRT:

Calreticulin

CYCP2;1:

Cyclin p2;1

DELTA-VPE:

Delta vacuolar processing enzyme

ER:

Endoplasmic reticulum

Ez.:

Enzyme

GO:

Gene Ontology

HS:

Heat stress

HSF:

Heat shock transcription factor

HSP:

Heat shock protein

HSR:

Heat-stress response

IGB:

Integrative Genome Browser

log2-CHG:

Log2-fold change

LTP:

Lipid transfer protein

Mbf1C:

Multiprotein Bridging Factor 1c

MIR:

MicroRNA

MS:

Murashige and Skoog

Mut:

cngc16

Norm:

Normal

PCA:

Principal component analysis

PE:

Paired-end

Q-PCR:

Quantitative PCR

ROS:

Reactive oxygen species

Sig. adj. p-val:

Significance based on adjusted p-value ≤0.01

SLP:

Calcineurin-like metallo-phosphoesterase superfamily protein

sncRNA:

Small non-coding RNA

SPL:

Squamosa Promoter Binding Protein-Like

T-DNA :

Transfer DNA

TF:

Transcription factor

UBQ7:

UBIQUITIN 7

UPR:

Unfolded protein response

VAMP:

Vesicle-associated membrane protein

WT:

Wild type

References

  1. Zinn KE, Tunc-Ozdemir M, Harper JF. Temperature stress and plant sexual reproduction: uncovering the weakest links. J Exp Bot. 2010;61:1959–68.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  2. Young LW, Wilen RW, Bonham-Smith PC. High temperature stress of Brassica napus during flowering reduces micro- and megagametophyte fertility, induces fruit abortion, and disrupts seed production. J of Exp Bot. 2004;55:485–95.

    Article  CAS  Google Scholar 

  3. Giorno F, Wolters-Arts M, Mariani C, Rieu I. Ensuring reproduction at high temperatures: the heat stress response during anther and pollen development. Plants. 2013;2:489–506.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  4. Hedhly A, Hormaza JI, Herrero M. Global warming and sexual plant reproduction. Trends Plant Sci. 2009;14:30–6.

    Article  PubMed  CAS  Google Scholar 

  5. Abiko M, Akibayashi K, Sakata T, Kimura M, Kihara M, Itoh K, Asamizu E, Sato S, Takahashi H, Higashitani A. High-temperature induction of male sterility during barley (Hordeum vulgare L.) anther development is mediated by transcriptional inhibition. Sex Plant Reprod. 2005;18:91–100.

    Article  CAS  Google Scholar 

  6. Kim SY, Hong CB, Lee I. Heat shock stress causes stage-specific male sterility in Arabidopsis thaliana. J Plant Res. 2001;114:301–7.

    Article  Google Scholar 

  7. Müller F, Rieu I. Acclimation to high temperature during pollen development. Plant Reprod. 2016;29:107–18.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  8. Fragkostefanakis S, Mesihovic A, Hu Y, Schleiff E. Unfolded protein response in pollen development and heat stress tolerance. Plant Reprod. 2016;29:81–91.

    Article  PubMed  CAS  Google Scholar 

  9. Deng Y, Srivastava R, Quilichini TD, Dong H, Bao Y, Horner HT, Howell SH. IRE1, a component of the unfolded protein response signaling pathway, protects pollen development in Arabidopsis from heat stress. Plant J. 2016;2:193–204.

    Article  CAS  Google Scholar 

  10. Mittler R, Finka A, Goloubinoff P. How do plants feel the heat? Trends Biochem Sci. 2012;37:118–25.

    Article  PubMed  CAS  Google Scholar 

  11. Ohama N, Sato H, Shinozaki K, Yamaguchi-Shinozaki K. Transcriptional regulatory network of plant heat stress response. Trends Plant Sci. 2017;22:53–65.

    Article  PubMed  CAS  Google Scholar 

  12. Rieu I, Twell D, Firon N. Pollen development at high temperature: from acclimation to collapse. Plant Physiol. 2017;173:1967–76.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Endo M, Tsuchiya T, Hamada K, Kawamura S, Yano K, Ohshima M, Higashitani A, Watanabe M, Kawagishi-Kobayashi M. High temperatures cause male sterility in rice plants with transcriptional alterations during pollen development. Plant Cell Physiol. 2009;50:1911–22.

    Article  PubMed  CAS  Google Scholar 

  14. Zhang X, Xiong H, Liu A, Zhou X, Peng Y, Li Z, Luo G, Tian Z, Chen X. Microarray data uncover the genome-wide gene expression patterns in response to heat stress in rice post-meiosis panicle. J Plant Biol. 2014;57:327–36.

    Article  CAS  Google Scholar 

  15. Yang J, Chen X, Zhu C, Peng X, He X, Fu J, Ouyang L, Bian J, Hu L, Sun X, Xu J, He H. RNA-Seq reveals differentially expressed genes of rice (Oryza sativa) spikelet in response to temperature interacting with nitrogen at meiosis stage. BMC Genomics. 2015;16:959–77.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Frank G, Pressman E, Ophir R, Althan L, Shaked R, Freedman M, Shen S, Firon N. Transcriptional profiling of maturing tomato (Solanum lycopersicum L.) microspores reveals the involvement of heat shock proteins, ROS scavengers, hormones, and sugars in the heat stress response. J Exp Bot. 2009;60:3891–908.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Bokszczanin KL, Krezdorn N, Fragkostefanakis S, Müller S, Rycak L, Chen Y, Hoffmeier K, Kreutz J, Paupière MJ, Chaturvedi P, Iannacone R, Müller F, Bostan H, Chiusano ML, Scharf K, Rotter B, Schleiff E, Winter P. Identification of novel small ncRNAs in pollen of tomato. BMC Genomics. 2015;16:714–33.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. Keller M, Hu Y, Mesihovic A, Fragkostefanakis S, Schleiff E, Simm S. Alternative splicing in tomato pollen in response to heat stress. Dna Res. 2017;24:205–17.

    Article  PubMed  CAS  Google Scholar 

  19. Loraine A, Blakley I, Jagadeesan S, Harper J, Miller G, Firon N. Analysis and visualization of RNA-Seq expression data using RStudio, Bioconductor, and integrated genome browser. Methods Mol Biol. 2015;1284:481–501.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Tunc-Ozdemir M, Tang C, Rahmati Ishka M, Brown E, Groves NR, Myers CT, Rato C, Poulsen LR, McDowell S, Miller G, Mittler R, Harper JF. A cyclic nucleotide-gated channel (CNGC16) in pollen is critical for stress tolerance in pollen reproductive development. Plant Physiol. 2013;16:1010–20.

    Article  CAS  Google Scholar 

  21. Cheng CY, Krishnakumar V, Chan AP, Thibaud-Nissen F, Schobel S, Town CD. Araport11: a complete reannotation of the Arabidopsis thaliana reference genome. Plant J. 2017;4:789–804.

    Article  CAS  Google Scholar 

  22. Loraine AN, McCormick S, Estrada A, Patel K, Qin P. RNA-Seq of Arabidopsis pollen uncovers novel transcription and alternative splicing. Plant Physiol. 2013;162:1092–109.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Qin Y, Leydon AR, Manziello A, Pandey R, Mount D, Denic S, Vasic B, Johnson MA, Palanivelu R. Penetration of the stigma and style elicits a novel transcriptome in pollen tubes, pointing to genes critical for growth in a pistil. PLoS Genetics. 2009;5:e1000621.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Tunc-Ozdemir M, Rato C, Brown E, Rogers S, Mooneyham A, Frietsch S, Myers CT, Poulsen LR, Malhó R, Harper JF. Cyclic nucleotide gated channels 7 and 8 are essential for male reproductive fertility. PLoS One. 2013;8:e55277.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Frietsch S, Wang YF, Sladek C, Poulsen LR, Romanowsky SM, Schroeder JI, Harper JF. A cyclic nucleotide-gated channel is essential for polarized tip growth of pollen. Proc Natl Acad Sci U S A. 2007;104:14531–6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Umate P. Genome-wide analysis of the family of light-harvesting chlorophyll a/b-binding proteins in Arabidopsis and rice. Plant Signal Behav. 2010;12:1537–42.

    Article  CAS  Google Scholar 

  27. Liu XL, Liu L, Niu QK, Xia C, Yang KZ, Li R, Chen LQ, Zhang XQ, Zhou Y, Ye D. Male gametophyte defective 4 encodes a rhamnogalacturonan II xylosyltransferase and is important for growth of pollen tubes and roots in Arabidopsis. Plant J. 2011;65:647–60.

    Article  PubMed  CAS  Google Scholar 

  28. Dai X, Zhao P. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Research. 2011; https://doi.org/10.1093/nar/GKR319.

  29. Stief A, Altmann S, Hoffmann K, Datt Pant B, Scheible W, Bäurle I. Arabidopsis miR156 regulates tolerance to recurring environmental stress through SPL transcription factors. Plant Cell. 2014;26:1792–807.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Li SB, Xie ZZ, Hu CG, Zhang JZ. A review of auxin response factors (ARFs) in plants. Front Plant Sci. 2016;7:47.

    PubMed  PubMed Central  Google Scholar 

  31. Zhang S, Liu Y, Yu B. New insights into pri-miRNA processing and accumulation in plants. Wiley Interdisciplinary Reviews: RNA. 2015;6:533–45.

    Article  PubMed  CAS  Google Scholar 

  32. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, Jensen LJ, von Mering C. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45:D362–8.

    Article  PubMed  CAS  Google Scholar 

  33. Schmid M, Davison TS, Henz SR, Pape UJ, Demar M, Vingron M, Schölkopf B, Weigel D, Lohmann JU. A gene expression map of Arabidopsis thaliana development. Nat Genet. 2005;37:501–6.

    Article  PubMed  CAS  Google Scholar 

  34. Wan XL, Zhou Q, Wang YY, Wang WN, Bao MZ, Zhang JW. Identification of heat-responsive genes in carnation (Dianthus caryophyllus L.) by RNA-Seq. Front Plant Sci. 2015;6:519–32.

    PubMed  PubMed Central  Google Scholar 

  35. Rizhsky L, Liang H, Mittler R. The combined effect of drought stress and heat shock on gene expression in tobacco. Plant Physiol. 2002;130:1143–51.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Rizhsky L, Liang HJ, Shuman J, Shulaev V, Davletova S, Mittler R. When defense pathways collide. The response of Arabidopsis to a combination of drought and heat stress. Plant Physiol. 2004;134:1683–96.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Bhardwaj AR, Joshi G, Kukreja B, Malik V, Arora P, Pandey R, Shukla RN, Bankar KG, Katiyar-Agarwal S, Goel S, Jagannath A, Kumar A, Agarwal M. Global insights into high temperature and drought stress regulated genes by RNA-Seq in economically important oilseed crop Brassica juncea. BMC Plant Biology. 2015;15:9–24.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Vihervaara A, Sergelius C, Vasara J, Blom M, Elsing A, Roos-Mattjus P, Sistonen L. Transcriptional response to stress in the dynamic chromatin environment of cycling and mitotic cells. Proc Natl Acad Sci U S A. 2013;110:E3388–97.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Nishizawa-Yokoi A, Yoshida E, Yabuta Y, Shigeoka S. Analysis of the regulation of target genes by an Arabidopsis heat shock transcription factor, HsfA2. Biosci Biotechnol Biochem. 2009;73:890–5.

    Article  PubMed  CAS  Google Scholar 

  40. Suzuki N, Bajad S, Shuman J, Shulaev V, Mittler R. The transcriptional co-activator MBF1c is a key regulator of thermotolerance in Arabidopsis thaliana. J Biol Chem. 2008;283:9269–75.

    Article  PubMed  CAS  Google Scholar 

  41. Kotak S, Larkindale J, Lee U, von Koskull-Döring P, Vierling E, Scharf KD. Complexity of the heat stress response in plants. Curr Opin Plant Biol. 2007;10:310–6.

    Article  PubMed  CAS  Google Scholar 

  42. Qu AL, Ding YF, Jiang Q, Zhu C. Molecular mechanisms of the plant heat stress response. Biochem Biophys Res Commun. 2013;432:203–7.

    Article  PubMed  CAS  Google Scholar 

  43. Zhong L, Zhou W, Wang H, Ding S, Lu Q, Wen X, Peng L, Zhang L, Lu C. Chloroplast small heat shock protein HSP21 interacts with plastid nucleoid protein pTAC5 and is essential for chloroplast development in Arabidopsis under heat stress. Plant Cell. 2013;25:2925–43.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Suzuki N, Mittler R. Reactive oxygen species and temperature stresses: a delicate balance between signaling and destruction. Physiologia Plantarum. 2006;126:45–518.

    Article  CAS  Google Scholar 

  45. Baxter A, Mittler R, Suzuki N. ROS as key players in plant stress signalling. J Exp Bot. 2014;65:1229–40.

    Article  PubMed  CAS  Google Scholar 

  46. Vanderauwera S, Suzuki N, Miller G, van de Cotte B, Morsa S, Ravanat JL, Hegie A, Triantaphylidès C, Shulaev V, Van Montagu MC, Van Breusegem F, Mittler R. Extranuclear protection of chromosomal DNA from oxidative stress. Proc Natl Acad Sci U S A. 2011;108:1711–6.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Swindell WR. The association among gene expression responses to nine abiotic stress treatments in Arabidopsis thaliana. Genetics. 2006;174:1811–24.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. Serivichyaswat P, Ryu HS, Kim W, Kim S, Chung KS, Kim JJ, Ahn JH. Expression of the floral repressor miRNA156 is positively regulated by the AGAMOUS-like proteins AGL15 and AGL18. Mol Cells. 2015;38:259–66.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Rowland O, Zheng H, Hepworth SR, Lam P, Jetter R, Kunst L. CER4 encodes an alcohol-forming fatty acyl-coenzyme a reductase involved in cuticular wax production in Arabidopsis. Plant Physiol. 2006;142:866–77.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Aarts MGM, Hodge R, Kalantidis K, Florack D, Wilson ZA, Mulligan BJ, Stiekema WJ. The Arabidopsis MALE STERILITY 2 protein shares similarity with reductases in elongation/condensation complexes. Plant J. 1997;12:615–23.

    Article  PubMed  CAS  Google Scholar 

  51. Fabrice TN, Vogler H, Draeger C, Munglani G, Gupta S, Herger AG, Knox P, Grossniklaus U, Ringli C. LRX proteins play a crucial role in pollen grain and pollen tube Cell Wall development. Plant Physiol. 2018;3:1981–92.

    Article  CAS  Google Scholar 

  52. Wang X, Wang K, Yin G, Liu X, Liu M, Cao N, Duan Y, Gao H, Wang W, Ge W, Wang J, Li R, Guo Y. Pollen-expressed leucine-rich repeat Extensins are essential for pollen germination and growth. Plant Physiol. 2018;3:1993–2006.

    Article  CAS  Google Scholar 

  53. Mecchia MA, Santos-Fernandez G, Duss NN, Somoza SC, Boisson-Dernier A, Valeria G, Martínez-Bernardini A, Fabrice TN, Ringli C, Muschietti JP, Grossniklaus U. RALF4/19 peptides interact with LRX proteins to control pollen tube growth in Arabidopsis. Science. 2017;358:1600–3.

    Article  PubMed  CAS  Google Scholar 

  54. Sede AR, Borassi C, Wengier DL, Mecchia MA, Estevez JM, Muschietti JP. Arabidopsis pollen extensins LRX are required for cell wall integrity during pollen tube growth. FEBS Letter. 2018;2:233–43.

    Article  CAS  Google Scholar 

  55. Gangadhar BH, Sajeesh K, Venkatesh J, Baskar V, Abhinandan K, Yu JW, Prasad R, Mishra RK. Enhanced tolerance of transgenic potato plants over-expressing non-specific lipid transfer Protein-1 (StnsLTP1) against multiple abiotic stresses. Front Plant Sci. 2016;7:1228.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Okazaki Y, Saito K. Roles of lipids as signaling molecules and mitigators during stress response in plants. Plant J. 2014;79:584–96.

    Article  PubMed  CAS  Google Scholar 

  57. McDowell SC, López-Marqués RL, Cohen T, Brown E, Rosenberg A, Palmgren MG, Harper JF. Loss of the Arabidopsis thaliana P4-ATPases ALA6 and ALA7 impairs pollen fitness and alters the pollen tube plasma membrane. Front Plant Sci. 2015;6:197.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Gao F, Han X, Wu J, Zheng S, Shang Z, Sun D, Zhou R, Li B. A heat-activated calcium-permeable channel-Arabidopsis cyclic nucleotide-gated ion channel 6- is involved in heat shock responses. Plant J. 2012;70:1056–69.

    Article  PubMed  CAS  Google Scholar 

  59. Finka A, Goloubinoff P. The CNGCb and CNGCd genes from Physcomitrella patens moss encode for thermosensory calcium channels responding to fluidity changes in the plasma membrane. Cell Stress Chaperones. 2014;19:83–90.

    Article  PubMed  CAS  Google Scholar 

  60. Charpentier M, Sun J, Vaz Martins T, Radhakrishnan GV, Findlay K, Soumpourou E, Thouin J, Véry AA, Sanders D, Morris RJ, Oldroyd GE. Nuclear-localized cyclic nucleotide-gated channels mediate symbiotic calcium oscillations. Science. 2016;352:1102–5.

    Article  PubMed  CAS  Google Scholar 

  61. Mi H, Muruganujan A, Casagrande JT, Thomas PD. Large-scale gene function analysis with the PANTHER classification system. Nat Protoc. 2013;8:1551–66.

    Article  PubMed  CAS  Google Scholar 

  62. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;15:2114–20.

    Article  CAS  Google Scholar 

  63. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  64. Swarbreck D, Wilks C, Lamesch P, Berardini TZ, Garcia-Hernandez M, Foerster H, Li D, Meyer T, Muller R, Ploetz L, Radenbaugh A, Singh S, Swing V, Tissier C, Zhang P, Huala E. The Arabidopsis information resource (TAIR): gene structure and function annotation. Nucleic Acids Res. 2008;36:1009–14.

    Article  CAS  Google Scholar 

  65. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;4:357–60.

    Article  CAS  Google Scholar 

  66. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;16:2078–9.

    Article  CAS  Google Scholar 

  67. Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;7:923–30.

    Article  CAS  Google Scholar 

  68. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15:550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  69. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B. 1995;57:289–300.

    Google Scholar 

  70. 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;4:402–8.

    Article  CAS  Google Scholar 

  71. Thimm O, Bläsing O, Gibon Y, Nagel A, Meyer S, Krüger P, Selbig J, Müller LA, Rhee SY, Stitt M. MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 2004;6:914–39.

    Article  CAS  Google Scholar 

  72. Borges F, Gomes G, Gardner R, Moreno N, McCormick S, Feijó JA, Becker JD. Comparative transcriptomics of Arabidopsis sperm cells. Plant Physiology. 2008;148:1168–81.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

Support in performing PANTHER GO analyses was provided by Drs. Won C. Yim and Juli Petereit.

Funding

This work was supported by grants to JFH and GM BARD IS-4652-13, and JFH from USDA Hatch and NSF IOS 1656774. The bioinformatics and statistics in this study were made possible by a grant from the National Institute of General Medical Sciences (GM103440) from the National Institutes of Health. The funders had no role in the designing the research, data collection, analysis, or manuscript preparation.

Availability of data and materials

The datasets generated during the current study are available with accession number SRP110833 in the NCBI Sequence Read Archive, https://www.ncbi.nlm.nih.gov/sra/.

Author information

Authors and Affiliations

  1. Department of Biochemistry and Molecular Biology, University of Nevada, Reno, MS330, Howard Building, Reno, NV, 89557, USA

    Maryam Rahmati Ishka, Elizabeth Brown, Chrystle Weigand, Karen A. Schlauch & Jeffrey F. Harper

  2. Nevada INBRE Bioinformatics Core, University of Nevada, Reno, Reno, NV, 89557, USA

    Richard L. Tillett & Karen A. Schlauch

  3. The Mina and Everard Goodman Faculty of Life Sciences Bar Ilan University, 52900, Ramat-Gan, Israel

    Gad Miller

Authors
  1. Maryam Rahmati Ishka
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Elizabeth Brown
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Chrystle Weigand
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Richard L. Tillett
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Karen A. Schlauch
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Gad Miller
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jeffrey F. Harper
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MRI designed and performed the experiments and wrote the manuscript; EB and CW assisted with experiments; RLT and KAS performed bioinformatics and statistical analyses; GM provided data analysis; JFH supervised the project. All authors read and approved the manuscript.

Corresponding author

Correspondence to Jeffrey F. Harper.

Ethics declarations

Ethics approval and consent to participate

All plant materials used were derived from Arabidopsis thaliana seed stocks available from the Arabidopsis Biological Resource Center at the Ohio State University (https://abrc.osu.edu). Research here did not involve field studies. All handling of plant materials complied with research guidelines at the University of Nevada, Reno.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

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

Additional files

Additional file 1:

Hot/Cold stress cycle. Diagram showing the Hot/Cold stress-cycle used here for growing plants from which pollen samples for RNA-Seq experiment were harvested at the end of HS-peak at 40 °C. See Methods and [20] for more details. (PPTX 50 kb)

Additional file 2:

Raw expression counts for WT and cngc16 pollen with and without HS. Read counts generated via FeatureCounts, see Methods) for all 12 replicates are shown before normalization and exclusions. (XLSX 4358 kb)

Additional file 3:

Normalized transcript expression counts for WT and cngc16 pollen with and without HS. Expression counts from Additional file 2 were subjected to exclusions/filters and normalized as described in methods. Expression data from other transcriptome studies (microarray and RNA-Seq) were added to the table for comparisons. These included an RNA-Seq data set for WT pollen and seedlings from Loraine et al. 2013 (PMCID: PMC3668042 [22]), two pollen microarray experiments from Qin et al. 2009 (PMCID: PMC2726614 [23]) and Borges et al. 2008 (PMCID: PMC2556834 [72]), and finally HS seedlings from Schmid et al. 2005 (PMID:15806101 [33]), respectively. Ratios of expression between pollen and seedling are based on Loraine et al. 2013 (PMCID: PMC3668042 [22]). In cases where the seedling value was below the limit of detection, a minimal value of 0.0019 was substituted in its place as a denominator (0.0019 was the RPKM for ATCG00860 and was the lowest value reported in Loraine et al. 2013 (PMCID: PMC3668042 [22]). Ratio of expression between semi-in vivo pollen tube over dry pollen is based on Qin et al. 2009 (PMCID: PMC2726614 [23]). HS dependent changes in transcript abundance in shoots were based on publicly available data using the AtGenExpress Visualization Tool (AVT) (http://jsp.weigelworld.org/expviz/expviz.jsp, Schmid et al. 2005 (PMID:15806101 [33] for seedlings exposed to one hour HS at 38 °C). The log2-fold change was calculated based on a comparison of means of normalized values for two heat-stressed and two non-stressed seedling samples. NA stands for not available. Not Calculated, refers to a value not being calculated because one of the input sample read counts was considered to have an extreme outlier (see AT2G42540 and ATMG01360). (XLSX 8952 kb)

Additional file 4:

Library size and principal component analysis. a. Table showing library sizes of each sample. b. A principal component analysis (PCA) of the filtered data showing that 87% of the variance of the samples can be explained by differences in the stress states. See methods for more details. Control and heat correspond to normal and HS conditions, respectively. (PPTX 43 kb)

Additional file 5:

A transcript profile comparison to evaluate purity of pollen samples used for RNA-Seq. A subset of 12 genes was used to compare relative purities of pollen samples in the current pollen transcriptome study to those from a RNA-Seq study from Loraine et al. [22] (yellow highlights) or a microarray experiment from Qin et al. 2009 [23] (purple highlights). Four references genes were chosen to generate normalization factors that could be used to adjust expression values in Loraine et al. [22] and Qin et al. 2009 [23] to allow a relative comparison of the three data sets for WT pollen under control (normal) conditions. For a control group, three CNGC genes were chosen that displayed low to moderate levels of expression (Tunc-Ozdemir et al. 2013 [24] and Frietsch et al. 2007 [25]). As markers for potential contamination from photosynthetic tissues, five different nuclear encoded genes were chosen that are associated with either photosystems I/II, or chlorophyll A-B binding proteins (Umate 2010 [26]). Average relative ratios are shown for each of the four different pollen samples in comparison to both Loraine et al. [22] and Qin et al. [23]. (XLSX 19 kb)

Additional file 6:

Integrated Genome Browser (IGB) screenshot showing cngc16 RNA-Seq reads primarily upstream of T-DNA insertion site. The green arrow identifies the position of T-DNA insertion in cngc16–2 (SAIL_726_B04). The observed reads aligning to cngc16 are primarily on the 5′ side of the T-DNA disruption site, with only a few reads observed at two disconnected downstream positions. This suggests that there were no detectable full-length transcripts. (PPTX 66 kb)

Additional file 7:

RNA-Seq validation using real-time Q-PCR. a. Comparison of expression values obtained from Q-PCR and RNA-Seq normalized to WT control (normal). The analysis was performed on two different reference genes separately (CYCP2 (AT3G21870) and UBQ7 (AT2G35635)). b. Primer sequences used for real-time Q-PCR. (XLSX 16 kb)

Additional file 8:

Predicted targets for HS-modulated microRNAs. Target predictions for microRNAs were conducted with psRNATarget (http://plantgrn.noble.org/psRNATarget/, Dia and Zhao 2011 (PMCID: PMC3125753 [28])) using Arabidopsis thaliana unigene library (2017 update). MIR156A and C have the same targets. ND stands for not detected in current pollen transcriptome study. NA stands for not applicable because no target was predicted using psTarget. (XLSX 57 kb)

Additional file 9:

Sorted lists of differentially expressed transcripts in WT and cngc16 mutant. HS-dependent changes were sorted and extracted from the entire transcriptome analysis in Additional file 3. For genes listed in a given category, corresponding transcripts showed ≥2-fold change with an adjusted p-value ≤0.01. a. Tally of 21 pre-existing transcriptome differences identified in a direct comparison between WT and cngc16 mutant pollen from unstressed (control/normal) plants. b. Tally of 471 HS-dependent transcriptome changes that were only classified as significant changes in WT. Genes shown include only WT-specific changes and do not include changes that were also observed to be in common with cngc16 (see Additional file 9d). c. Tally of 2305 HS-dependent transcriptome changes that were only classified as significant changes in cngc16. Genes shown include only cngc16-specific changes and do not include changes that were also observed to be in common with WT (see Additional file 9d). d. Tally of 1631 HS-dependent transcriptome changes observed in both WT and cngc16. Genes shown include only common changes and do not include changes that were categorized as specific to WT or cngc16 (see Additional file 9b and c, respectively). e. 192 transcriptome differences from a direct comparison between heat-stressed WT and cngc16 mutant pollen. (XLSX 1638 kb)

Additional file 10:

Association network for HS-dependent changes in WT that were most different in cngc16 under HS. The figure was generated using STRING [32] and the top 23 genes identified in Table 3. (XLSX 154 kb)

Additional file 11:

HS-dependent transcript abundance changes corresponding to transcription factors in WT and cngc16 pollen. Differentially expressed transcription factors (TFs) are organized to show HS-dependent differences that are WT-specific (23 changes, yellow highlight), cngc16-specific (98 changes, blue highlight), or common to both mutant and WT (66 changes, white highlight). a A comparison is provided to HS-dependent changes observed for aerial parts of seedlings exposed to a one hour 38 °C HS based on publicly available microarray data in AtGenExpress, (http://jsp.weigelworld.org/expviz/expviz.jsp; Schmid et al. 2005 (PMID: 15806101 [34]). The log2-fold change was calculated based on means of normalized values for two heat-stressed seedlings compared to two non-stressed seedlings. NA stands for not available because of absence of probeset on microarray experiment. (XLSX 63 kb)

Additional file 12:

230 Genes with putative Ca2+-signaling related functions. 230 genes with putative Ca2+-signaling related functions were identified here with quantifiable expression levels in the pollen RNA-Seq data set. Of those 230, 40 genes (yellow highlights) showed HS-dependent changes (≥ 2-fold changes with adjusted p-value ≤0.01) in WT pollen. Similar HS-dependent changes were also observed for cngc16 mutant, with seven potential exceptions with lower degrees of significance shown in bold. a Ratio of expression between pollen and seedling is based on Loraine et al. 2013 (PMCID: PMC3668042 [22]). b A comparison is provided to HS-dependent changes observed for aerial parts of seedlings exposed to a one hour 38 °C HS based on publicly available microarray data in AtGenExpress, (http://jsp.weigelworld.org/expviz/expviz.jsp; Schmid et al. 2005 (PMID: 15806101 [33]). The log2-fold change was calculated based on means of normalized values for two heat-stressed seedlings compared to two non-stressed seedlings. (XLSX 15 kb)

Additional file 13:

A comparison of HS-dependent changes in pollen to 67 multi-stress response genes in vegetative tissues. From a list of 67 multi-stress genes curated by Swindell 2006 (PMCID: PMC1698639 [47]; highlighted in purple), 19 genes showed detectable expression in pollen. Among those, only three genes showed significant changes in pollen HS (red font). (XLSX 596 kb)

Additional file 14:

GO analyses on HS-dependent changes in WT and cngc16. a. The number of genes in each GO category is shown for genes with HS-dependent changes observed in both WT (green header) and cngc16 (orange header). Enrichment above an expected is shown along with a p-value. In a simple contrast analysis, a ratio of gene numbers in cngc16 and WT is calculated for each GO category. The analysis was done as a PANTHER Overrepresentation Test (release 2017–04-13 [61]) using a GO Ontology database (released 2017–08-14) with 27,060 reference genes for Arabidopsis thaliana. NA stands for not applicable because no genes were detected in this category for either WT or cngc16. b. Uploads used for HS-dependent changes with ≥2-fold changes and adjusted p-value ≤0.01. (PPTX 138 kb)

Additional file 15:

GO analysis on the 192 largest differences between WT and cngc16 under HS. A GO analysis pie chart is shown for Molecular Function (a), Cellular Component (b), and Biological Process (c) generated using an upload of Additional file 3 or Additional file 9e column R listing the differences (≥ 2-fold and adjusted p-value ≤0.01) between WT and cngc16 HS-transcriptomes. Categories were defined using PANTHER Overrepresentation Test (release 2017–04-13 [61]) using a GO Ontology database (released 2017–08-14) with 27,060 reference genes for Arabidopsis thaliana. Gene categories shown displayed enrichments with a p-value of ≤0.05. (PPTX 993 kb)

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

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rahmati Ishka, M., Brown, E., Weigand, C. et al. A comparison of heat-stress transcriptome changes between wild-type Arabidopsis pollen and a heat-sensitive mutant harboring a knockout of cyclic nucleotide-gated cation channel 16 (cngc16). BMC Genomics 19, 549 (2018). https://doi.org/10.1186/s12864-018-4930-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-4930-4

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us