Codominant grasses differ in gene expression under experimental climate extremes in native tallgrass prairie

Site description and experimental treatments

The study was carried out within the context of an existing long-term climate change experiment, the Rainfall Manipulation Plots (RaMPs), located at the Konza Prairie Biological Station in north-eastern Kansas (39°05′N, 96°35′W). Kansas State University, Manhattan, KS, USA and the Konza Prairie Biological Station granted explicit permission to the authors to sample with minimal impact within the RaMPs. The RaMPs is located in a native, annually burned site and consists of twelve 14 × 9 m greenhouse shelters (without walls) equipped with a clear (UV transparent) polyethylene roof to exclude natural rainfall inputs (Fay et al., 2011). Our experimental plots were located in two RaMPs (RaMP 12 and 13) in areas outside the 6 × 6 m experimental plots, but still located underneath the shelter infrastructure. Each of these areas is approximately 3 × 8 m in size, within which we located a 3 × 6 m experimental sampling plot. The RaMP 12 sampling plot was watered from late-May to mid-Aug to create a watered condition, whereas all ambient rainfall was excluded from the RaMP 13 sampling plot to create a drought. For both the watered and drought plots, a controlled high heat treatment was achieved by installing pairs of rectangular infrared heating lamps (Kalglo 2000 W; Kalglo Electronic Co Inc., Bethlehem, PA, USA) (Fig. S1). This resulted in a high heat treatment zone with a daytime target maximum of +8 °C above ambient midday temperature (Fig. S2), alongside ambient temperature treatment zones. The four treatments allowed us to examine the effects of drought and heat individually along with their interaction. The high heat treatment was imposed for an 18-day period (July 17 to August 4), when heat waves have generally occurred in the past (Hoover, Knapp & Smith, 2014b).

Prior to initiation of the experiment, canopy temperature in the watered sampling plot was measured using an infrared thermometer mounted on a movable platform (approx. 0.5 m above the canopy). Soil moisture was monitored at a depth of 0–15 cm with 30-cm time-domain reflectometery probes (Model CS616, Campbell Scientific, Logan, Utah, USA) inserted at a 45° angle (see Supplemental Information).

Plant sampling and measurements

The focal species, A. gerardii and S. nutans, are both are rhizomatous C₄ grasses that reproduce primarily vegetatively via belowground buds on rhizomes (Brejda, Moser & Waller, 1989; Carter & VanderWeide, 2014), which form dense intermixed stands, making it virtually impossible to differentiate between clones in the field (Avolio, Chang & Smith, 2011). We sampled individuals of A. gerardii and S. nutans from native populations growing within the experimental treatment plots during two sampling campaigns conducted at Day 4 and Day 18 of the heat wave. Each sampling campaign was conducted between 11:00 and 15:00 CDT to allow for collection of leaf temperature and water status (see below).

During each sampling campaign, we sampled two, morphological similar individuals (tiller or ramets, with 3–5 fully expanded leaves) of each species within the high heat zone and ambient temperature zone in both the watered and drought sampling plots (n = 2 samples per species, four treatments, and two campaign dates, or n = 16 per species, N = 32 total samples). While a sample size of two per species and treatment combination is relatively small, we believe this sample size was appropriate given that our focus was on broadly detecting interspecific differences under the high heat and drought conditions. Although we did not control for plant genotype, we collected our samples within a limited sampling area (10 × 10 cm) to minimize genotypic differences among samples. Leaf tissue was collected from individuals located within each treatment within a five-minute window. For each individual, the first or second fully expanded leaf was randomly selected for genomic analysis to ensure similar leaf age. The entire leaf was clipped and immediately flash-frozen and stored in liquid nitrogen until brought to the laboratory. Immediately after, we measured leaf temperature (T_leaf) and midday leaf water potential (Ψ_mid) on the remaining fully expanded leaf. T_leaf was measured using a LI-6400 system (LiCOR, Inc., Lincoln, NE, USA). The whole leaf was then collected for determination of mid-day leaf water potential (LWP) using an Scholander-type pressure chamber (PMS Instruments, Inc., Corvallis, OR, USA).

RNA preparation and microarray hybridization

Leaf tissue samples were stored in an −80 freezer prior to RNA extraction. Total RNA was extracted from the 32 leaf samples for both species using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) (McCarty, 1986), and further purified with the RNeasy kit (Invitrogen, Carlsbad, CA, USA). RNA quantity was measured by a NanoDrop spectrophotometer (Nanodrop products, Thermo Scientific, Wilmington, DE, USA). The verification of RNA quality, preparation of cDNA, and the subsequent steps leading to hybridization and array scanning were performed by Biotechnology Resources of Keck facility at Yale University (http://keck.med.yale.edu/). We used maize spotted cDNA arrays (SAM 1.2, GEO platform GPL4521) produced by the Center for Plant Genomics at Iowa State University for hybridization. The arrays included 15,680 maize cDNA probes (14,118 informative) isolated from maize ear tissue.

Quality control of heterologous hybridizations

In total, there were eight hybridizations for each species per sampling campaign (Table S1). Array image data were collected using GenePix software (Version 6; Axon, Downingtown, PA, USA). Prior to normalization across arrays, features with obvious abnormality and saturated signal were flagged and excluded from statistical analysis. Two steps were taken to minimize the probability of mishybridization and sequence divergence between the focal species and the model species (Leakey et al., 2009). First, we used stringent criteria by excluding spots with signal to noise ratios less than 3 or larger than 10 to decrease the inclusion of cross-hybridization artefacts (Verdnik, Handran & Pickett, 2002). Second, the cDNA sequences of the maize microarray SAM1.2 (18,862 sequences) were aligned against the de novo RNA-seq transcriptome data sets of A. gerardii and S. nutans (Hoffman & Smith, 2017), previously generated using Trinity (version 2.1.1, Haas et al., 2013). We only included BLASTN (Altschul et al., 1990; Altschul et al., 1997) hits with an e-value cutoff of 1e⁻¹⁰ and alignment length larger than 150 base pairs from the A. gerardii and S. nutans transcript data sets. After these two steps, 7,964 and 6,035 probe sequences were included in the analysis, accounting for 61.4% and 56.6% of the maize SAM 1.2 array probes for A. gerardii and S. nutans respectively. A total of 5,109 features were common to both species. Because features were screened by both the intensity of hybridization signal and sequence similarity, the intensity values of the included features were reliable for further expression analysis. These same techniques have also been validated previously using quantitative real-time PCR (qPCR) (Smith, Hoffman & Avolio, 2016).

Array data normalization and statistical analysis

An important source of systematic errors in two-color microarray experiments is the different properties of the dyes used to label the two samples (Tseng et al., 2001; Yang et al., 2001; Yang et al., 2002) and the hybridization variability from array to array. We used dye-swap design for the same pair of samples in the hybridizations (Table S1) to account for the dye effect (Dabney & Storey, 2007). Background signals were removed from median signal intensity and modelled similarly to Travers et al. (2010) to remove the array and dye effect: ${\displaystyle y_{ijk}=A_i+D_j+A_i{D}_j+{\epsilon }_{ijk},}$ where y is the median intensity for the k th gene on each array (i) with each dye (j), A is the array effect for each array (i), D is the dye effect for each dye (j), AD is the array × dye interaction, and ε_ijk is the stochastic error. Residuals from this model were adjusted by the minimum value to produce all positive residuals. To examine overall statistical effects, we used the residuals in the following model: ${\displaystyle r_{klmno}=S_l+W_m+T_n+C_o+S_l{W}_m+S_l{T}_n+W_m{T}_n+{\epsilon }_{klmno},}$ where r is the residual for each gene (k) with each species (l), water treatment (m), temperature (n), and sampling date (o), S is the species effect, W is the water treatment effect (plot), T is the temperature effect, and C is the sampling date effect. Residuals were used to generate log₂ expression ratios for the four variables: species (A. gerardii/S. nutans), water treatment (watered/droughted), temperature (ambient/heated), date (day 4/day 18). Any genes with missing signals were removed. We plotted the log₂ expression ratio against the log₁₀ intensity for each gene and performed a loess correction to normalize each set of log₂ values (Fig. S3). Then, for each gene without missing values, a linear model was performed to test each main effect (species, water treatment, temperature, and date) as well as selected interactions (species × water treatment, species × temperature, and water treatment × temperature). Because of the variation in genes present across arrays, each model was constructed only if appropriate data was present. In other words, to test species effect, both species had to express the given gene. P-values were adjusted using a Bonferroni correction to account false discovery across multiple tests. All analyses were performed using R (version 3.3.2).

Functional annotation, enrichment, and clustering

The functional annotation of transcripts was based on the Trinotate pipeline (version 3.0.1). We matched microarray probe sequences to known sequences using BLAST against the SwissProt annotated database (Apweiler et al., 2004), identified protein sequence homology using HMMER and Pfam (Finn, Clements & Eddy, 2011; Finn et al., 2015), and searched for known annotations within eggNOG and GO databases (The Gene Ontology Consortium, 2015; Huerta-Cepas et al., 2016). Ontology enrichment was determined using GOSeq (version 3.4, Young et al., 2010), a statistical package for R which accounts for multiple testing as well as differing probe lengths. Finally, clustering of gene modules was performed using the WGCNA package for R (version 1.51, Langfelder & Horvath, 2008) with a minimum module size of five genes.

Efficacy of the heat wave and drought treatments and impacts on T_leaf and Ψ_mid

On average, the heated (heat wave) treatment resulted in an 8 °C increase in canopy temperature (Fig. S2A). On average, the drought treatment decreased volumetric soil water content from 28% to 24% midway through the heat wave (day 9). The high heat treatment further decreased soil water content by 2% for the watered and 5% for the drought treatments (Fig. S2B). Overall, the combined effect of drought and heat resulted in a drop from 29% to 22% volumetric water content. The increase in canopy temperature with the high heat treatment was reflected in greater leaf temperature (T_leaf) for both species; A. gerardii and S. nutans had significantly higher T_leaf at both day 4 and 18 of the heat wave (Fig. S4). water content with the drought and high heat treatment were reflected in greater water stress in both species (i.e., more negative Ψ_mid, Fig. S4). For A. gerardii, the high heat treatment caused a large decrease in Ψ_mid, with this decline greatest at day 4 of the heat wave combined with drought (−0.9 MPa, Fig. S4). The decrease in Ψ_mid with the high heat treatment was most pronounced in S. nutans after 18 days of heat wave under drought (−1.7 MPa, Fig. S4).

Environment affects gene regulation in A. gerardii and S. nutans

Overall, 1,131 genes were shared across both species, 1,515 were shared across water treatment, 1,653 were shared across temperature treatment, and 1,390 were shared across date. Species (p < 0.001), water treatment (p < 0.001), and their interaction (p < 0.001) most significantly impacted gene expression. In other words, species gene expression response strongly depended on the drought environment. Temperature was only a weakly significant predictor of gene expression (p = 0.048) with no significant species by temperature interaction. Gene expression did not vary across sample date/duration of the heat wave.

Overall differences between A. gerardii and S. nutans

Of 1,131 genes found in both species, 160 differed significantly in their regulation between species. Genes with greater expression in A. gerardii were enriched in cellular metabolic process, biological regulation, and protein metabolic process, while genes with greater expression in S. nutans were enriched in negative regulation of metabolism, biological, and cellular processes, macromolecule catabolic process, and protein kinase activity (Fig. 1). Within cellular metabolic process, the most extreme differences were found in a methyltransferase and other transferases, GTP binding protein, Dihydrouridine synthase (Dus), as well as several transcription factors (Table 1). Among biological regulation genes, several transcription factors were strongly upregulated in A. gerardii. Protein metabolic processes included several ribosomal-related genes as well as fibrillarin upregulated in A. gerardii. Within genes significantly upregulated in S. nutans, the negative regulation (inhibition) category consisted of a finger protein as well as several membrane proteins like CMP-sialic acid transporter homolog (Table 1). Macromolecule catabolism genes included several proteasomes, 1,2-alpha-mannosidase, and a ubiquitin-conjugating enzyme. Among genes annotating to the term “stress”, 18 were upregulated in S. nutans versus 31 upregulated in A. gerardii. Genes annotating broadly to “regulation” showed 91 upregulated in A. gerardii versus 74 in S. nutans.

Gene clustering was performed for day 18 samples to detect species differences for both plots at the end of the heat wave. Similarly regulated modules or groups of genes may lead to a greater understanding of gene networks contributing to different species responses. One gene module significantly explained species differences in the watered treatment (p < 0.001, Fig. 2A) with genes generally expressed more highly in S. nutans. Two gene modules significantly explained species differences in the drought treatment (p = 0.01, Fig. 2B and p = 0.02, Fig. 2C respectively). Under drought, genes generally had lower expression in S. nutans.

Genes regulated in A. gerardii

In A. gerardii, 61 genes were significantly regulated in response to drought (5% of 1,148 total genes), with 24 genes upregulated under watered conditions and 37 upregulated under drought conditions. Few GO categories had strong enrichment (i.e., few genes per category). The drought treatment showed enrichment in response to osmotic stress, chromatin silencing, and lysosome. The watered treatment suggested greater abundance of xylose metabolism, sucrose metabolism, and ion transport (although each group contained only one gene) (Fig. 3A). Osmotic stress genes included an RNA-binding protein, ribosomal protein S3, and aconitate hydratase (Table 1). Within chromatin silencing genes, two histone acetyltransferases were upregulated under drought. Among all genes, 24 genes annotating to “stress” were upregulated in the watered treatment, versus 29 under drought. Only two genes (both within A. gerardii) responded significantly to temperature. One gene was upregulated in response to higher temperatures (Hsp70 protein); another was downregulated under higher temperatures (high mobility group-box domain).

Genes regulated in S. nutans

Sorghastrum nutans regulated more genes in response to drought than A. gerardii (23% of 762 genes total). Of these, 92 showed greater expression in the watered treatment while 82 showed greater expression under drought. Genes upregulated in the watered treatment showed GO enrichment in non-coding RNA (ncRNA) and RNA metabolism and nitrogen response. Genes upregulated under drought showed enrichment in response to water stress, external encapsulating structure, organophosphate metabolism, and cellular catabolism (Fig. 3B). Within the watered treatment ncRNA metabolism genes including ERBB-3 binding ribonuleoprotein, serrate RNA effector molecule, and pseudouridine synthase were upregulated (Table 1). Sorghastrum nutans in the watered treatment also showed greater expression of aquaporin NIP3-1, NEP1-interacting protein, and a transcriptional corepressor.

In contrast, S. nutans under drought showed greater expression of osmotic stress genes E3 ubiquitin ligase SUD1, 9 aldo-keto reductase, and hydrophobic protein LTI6A (Table 1). Among encapsulating structures, CMP-sialic acid transporter homolog, phosphatidylinositol kinase, pectin acetylesterase 8, and two glucuronosyltransferases (ranged from fold change of −1.47 to −1.71) were upregulated under drought. Catabolism genes within the drought treatment included 26S protease, DNA-directed RNA polymerase II Rpb7p, and phosphatidylinositol kinase. Lastly, the drought treatment showed increased expression of organophosphate metabolism genes including GDP-mannose 4,6 dehydratase, triosephosphate isomerase and phosphatidylinositol-4-phosphate 5-kinase. Among all genes, 12 (1.5%) genes annotating to “stress” were upregulated in the watered treatment, versus 20 (2.6%) under drought.