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)

In flowering plants, the male gametophyte (pollen) is one of the most vulnerable cells to temperature stress. In Arabidopsis thaliana, a pollen-specific CyclicNucleotide-Gated cationChannel 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.


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.
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_nonco-ding_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.
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   HS-dependent changes in transcript abundance in seedlings (aerial parts) based on publically available data using the AtGenExpress Visualization Tool (AVT) (http://jsp.weigelworld.org/expviz/expviz.jsp) [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 c Ratio of expression between pollen and seedling is based on [22]. NA stands for not applicable because it was not possible to calculate based on the information available d Because seedling value was below the limit of detection (< 0.0019), a minimal value of 0.0019 was used in its place as a denominator (0.0019 was the RPKM for ATCG00860 and was the lowest value reported in Loraine [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).

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).
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 Network group analysis was done using STRING [32] b HS-dependent changes in transcript abundance in seedlings (aerial parts) based on publically available data using the AtGenExpress Visualization Tool (AVT) (http://jsp.weigelworld.org/expviz/expviz.jsp) [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 c Ratio of expression between pollen and seedling is based on [22]. NA stands for not applicable because it was not possible to calculate based on the information available Functional annotation based on S (STRING, [32]), MM (MapMan, [71]), and T (TAIR, [64]), as noted 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 Ca 2+ -signaling related genes
Because CNGC16 is thought to function in both creating and sensing Ca 2+ -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 Ca 2+ -signaling (Ca 2+ -signaling related genes). Pollen expression was detected for 230 Ca 2+ -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 Ca 2+ -signaling related genes, the loss of CNGC16 did not have a significant impact. To evaluate whether Ca 2+ -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 Ca 2+ -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 Ca 2+ -related genes detected in pollen and seedlings.
Log2-Fold change in response to heat stress

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].
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 Ca 2+ -signals (Additional file 12). Of 40 Ca 2+ -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 Ca 2+ signal. While CNGC16's ion selectivity properties have not been corroborated, a Ca 2+ conductance has been suggested based on analogy to other closely related homologs [20]. For example, there is electrophysiology evidence for Ca 2+ conductance from a study on a closely related cngc6 knockout [58]. In this case, CNGC6-mediated Ca 2+ -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 Ca 2+ 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 Ca 2+ oscillation associated with the nucleus [60]. These Medicago homologs appear to be localized to the nuclear envelope, raising the potential that a HS-dependent Ca 2+ 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 Ca 2+ 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.
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.

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 Feature-Counts 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_non-coding_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].

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) 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