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)
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.
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.
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.
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) . 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.  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 ). 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.
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)  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  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 : 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 . 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  or ), 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.  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 . 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.
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 . 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 ). 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/, ) (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  and regulate floral development . 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 . 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).
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 , 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).
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.
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; ). 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.
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  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; ).
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 . 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 , and is essential for pollen fertility , 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 . In addition, there are many studies showing the importance of regulating lipid biogenesis pathways for cells to adapt to temperature stresses . As a specific example for pollen, a double knockout of lipid fippases ala6 and ala7 in Arabidopsis results in a stress-dependent pollen sterility  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 . For example, there is electrophysiology evidence for Ca2+ conductance from a study on a closely related cngc6 knockout . 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 . 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 . 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 ) (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.
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.
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.
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 . 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.
The stress regime with hot days and cold nights was performed as described in  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.
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  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  with BLASTN search via the command-line ncbi-blast+ (v2.5.0) application . 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  of the TAIR10 Arabidopsis thaliana reference genome sequence , 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) , using the spliced-read alignment tool HISAT2, version 2.0.5 . Alignments produced by HISAT2 were converted to the binary BAM format and sorted with samtools, version 1.3.1 . 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 . 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 . 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 . 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 . 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  (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 . 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  and MapMan version 3.6.0..
Auxin Response Factor
BCL-2-associated athanogene 6
- Ca2+ :
Calcineurin B subunit-like protein
Ca2+-binding endonuclease/exonuclease/phosphatase family
Calcineurin B-like (CBL)-interacting protein kinase
Cyclic nucleotide-gated cation channel
- Col-0 WT:
Columbia-0 wild type
Ca2+-dependent protein kinases
Delta vacuolar processing enzyme
Heat shock transcription factor
Heat shock protein
Integrative Genome Browser
Lipid transfer protein
Multiprotein Bridging Factor 1c
Murashige and Skoog
Principal component analysis
Reactive oxygen species
- Sig. adj. p-val:
Significance based on adjusted p-value ≤0.01
Calcineurin-like metallo-phosphoesterase superfamily protein
Small non-coding RNA
Squamosa Promoter Binding Protein-Like
- T-DNA :
Unfolded protein response
Vesicle-associated membrane protein
Zinn KE, Tunc-Ozdemir M, Harper JF. Temperature stress and plant sexual reproduction: uncovering the weakest links. J Exp Bot. 2010;61:1959–68.
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.
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.
Hedhly A, Hormaza JI, Herrero M. Global warming and sexual plant reproduction. Trends Plant Sci. 2009;14:30–6.
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.
Kim SY, Hong CB, Lee I. Heat shock stress causes stage-specific male sterility in Arabidopsis thaliana. J Plant Res. 2001;114:301–7.
Müller F, Rieu I. Acclimation to high temperature during pollen development. Plant Reprod. 2016;29:107–18.
Fragkostefanakis S, Mesihovic A, Hu Y, Schleiff E. Unfolded protein response in pollen development and heat stress tolerance. Plant Reprod. 2016;29:81–91.
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.
Mittler R, Finka A, Goloubinoff P. How do plants feel the heat? Trends Biochem Sci. 2012;37:118–25.
Ohama N, Sato H, Shinozaki K, Yamaguchi-Shinozaki K. Transcriptional regulatory network of plant heat stress response. Trends Plant Sci. 2017;22:53–65.
Rieu I, Twell D, Firon N. Pollen development at high temperature: from acclimation to collapse. Plant Physiol. 2017;173:1967–76.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Dai X, Zhao P. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Research. 2011; https://doi.org/10.1093/nar/GKR319.
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.
Li SB, Xie ZZ, Hu CG, Zhang JZ. A review of auxin response factors (ARFs) in plants. Front Plant Sci. 2016;7:47.
Zhang S, Liu Y, Yu B. New insights into pri-miRNA processing and accumulation in plants. Wiley Interdisciplinary Reviews: RNA. 2015;6:533–45.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Qu AL, Ding YF, Jiang Q, Zhu C. Molecular mechanisms of the plant heat stress response. Biochem Biophys Res Commun. 2013;432:203–7.
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.
Suzuki N, Mittler R. Reactive oxygen species and temperature stresses: a delicate balance between signaling and destruction. Physiologia Plantarum. 2006;126:45–518.
Baxter A, Mittler R, Suzuki N. ROS as key players in plant stress signalling. J Exp Bot. 2014;65:1229–40.
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.
Swindell WR. The association among gene expression responses to nine abiotic stress treatments in Arabidopsis thaliana. Genetics. 2006;174:1811–24.
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.
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.
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.
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.
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.
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.
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.
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.
Okazaki Y, Saito K. Roles of lipids as signaling molecules and mitigators during stress response in plants. Plant J. 2014;79:584–96.
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.
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.
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.
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.
Mi H, Muruganujan A, Casagrande JT, Thomas PD. Large-scale gene function analysis with the PANTHER classification system. Nat Protoc. 2013;8:1551–66.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;15:2114–20.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.
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.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;4:357–60.
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.
Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;7:923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15:550.
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.
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.
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.
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.
Support in performing PANTHER GO analyses was provided by Drs. Won C. Yim and Juli Petereit.
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/.
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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)
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 ), two pollen microarray experiments from Qin et al. 2009 (PMCID: PMC2726614 ) and Borges et al. 2008 (PMCID: PMC2556834 ), and finally HS seedlings from Schmid et al. 2005 (PMID:15806101 ), respectively. Ratios of expression between pollen and seedling are based on Loraine et al. 2013 (PMCID: PMC3668042 ). 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 ). Ratio of expression between semi-in vivo pollen tube over dry pollen is based on Qin et al. 2009 (PMCID: PMC2726614 ). 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  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)
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)
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.  (yellow highlights) or a microarray experiment from Qin et al. 2009  (purple highlights). Four references genes were chosen to generate normalization factors that could be used to adjust expression values in Loraine et al.  and Qin et al. 2009  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  and Frietsch et al. 2007 ). 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 ). Average relative ratios are shown for each of the four different pollen samples in comparison to both Loraine et al.  and Qin et al. . (XLSX 19 kb)
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)
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)
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 )) 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)
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)
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 ). 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)
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 ). 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 ). 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)
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 ; 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)
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 ) 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)
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 ) 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)
About this article
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