Genome-Wide Expression and Alternative Splicing in Domesticated Sunflowers (Helianthus annuus L.) under Flooding Stress

Approximately 10% of agricultural land is subject to periodic flooding, which reduces the growth, survivorship, and yield of most crops, reinforcing the need to understand and enhance flooding resistance in our crops. Here, we generated RNA-Seq data from leaf and root tissue of domesticated sunflower to explore differences in gene expression and alternative splicing (AS) between a resistant and susceptible cultivar under both flooding and control conditions and at three time points. Using a combination of mixed model and gene co-expression analyses, we were able to separate general responses of sunflower to flooding stress from those that contribute to the greater tolerance of the resistant line. Both cultivars responded to flooding stress by upregulating expression levels of known submergence responsive genes, such as alcohol dehydrogenases, and slowing metabolism-related activities. Differential AS reinforced expression differences, with reduced AS frequencies typically observed for genes with upregulated expression. Significant differences were found between the genotypes, including earlier and stronger upregulation of the alcohol fermentation pathway and a more rapid return to pre-flooding gene expression levels in the resistant genotype. Our results show how changes in the timing of gene expression following both the induction of flooding and release from flooding stress contribute to increased flooding tolerance.


Introduction
Historically, people have inhabited regions with abundant water for multiple reasons, but most of all, for successful agriculture. However, the association with water means that flooding typically is common as well, which can result in significant crop losses, especially for crops that lack resistance to flooding. Model-based estimates of crop losses resulting from flooding stress suggest that such losses are likely to increase dramatically with climate change due to the increased frequency and magnitude of coastal storms [1]. Moreover, there is some evidence that such predicted impacts are already being observed [2]. Thus, there is a need to better understand how some plants are able to resist flooding, with the long-term goal of developing crops with greater resilience to flooding stress.
From a plant physiological standpoint, flooding results in reduced gas exchange in submerged roots and shoots, which can slow down rates of respiration and photosynthesis [3] or lead to damage at the cellular level (e.g., in lipids, proteins and DNA) due to increased levels of reactive oxygen species (ROS) [4]. Plants that are frequently exposed to flooding, however, have developed several strategies to respond and adapt to flooding stress, both physiologically and morphologically. These strategies include quiescence, anaerobic energy production, rapid shoot growth to escape flooding, and the development of aerenchyma to facilitate internal gas exchange.
In the quiescence or "sit and wait" strategy, plants conserve energy such as ATP and carbohydrates to survive while submerged and to recover growth once flooding ends [3,5]. For example, tolerance to water submergence was shown to be negatively correlated with growth in 86 Arabidopsis accessions with natural variation for flooding tolerance, consistent with the quiescence strategy [6].
When oxygen, the essential substrate for respiration, becomes limited, plants are no longer able to produce enough ATPs to meet cellular demands [7]. Under oxygen limiting conditions, regeneration of nicotinamide adenine dinucleotide (NAD) occurs via two fermentation pathways: (i) lactic acid fermentation and ii) alcohol fermentation [7][8][9]. Lactic acid fermentation results in the potentially damaging accumulation of lactate, whereas alcohol (ethanol) fermentation is considered to be less harmful. Therefore, the extent of ethanol production is thought to be positively correlated with a plant's tolerance to flooding stress [10].
Wide-ranging development of air-filled aerenchyma tissue can also contribute to flooding resistance by facilitating gas exchange between the shoot and the root [11]. However, this mechanism is ineffective for complete submergence, which typically triggers rapid shoot growth to reach contact with air. The latter is referred to as an "escape strategy" [12] and occurs in both wild and crop plants [3,13,14]. The development of adventitious roots containing aerenchyma represents another escape strategy by minimizing the distance between tissues that are in contact with air (above the waterline) and root tips [14][15][16]. Adventitious root growth is promoted by increased ethylene gas [17,18] and is mediated by a combination of hormonal signals and reactive oxygen species [19]. The timing of adventitious root emergence is determined by the degree of flooding stress and the development stage of the plant, and it varies across plant species [16].
Transcriptomic analysis can provide additional clues regarding the mechanisms underlying resistance to flooding stress [20,21]. An example comes from soybean, which displays extensive variation for waterlogging tolerance [22][23][24][25][26]. Transcriptional responses were explored by RNA-Seq, and flooding stress was associated with downregulation of photosynthesis and chlorophyll synthesis related genes, consistent with reduced photosynthetic efficiency (i.e., quiescence) observed at the whole plant level [27]. Likewise, a combination of transcriptome and proteome analyses showed that alcohol fermentation related genes were significantly upregulated in soybean under flooding stress [28], indicating a role for fermentation as well. Arora et al. [29] conducted RNA-Seq in a flooding tolerant maize inbred line under waterlogging and identified differentially expressed genes involved in several key pathways, including energy-production, programmed cell death, aerenchyma formation, and ethylene responsiveness, suggesting that multiple mechanisms underlie responses to flooding stress. There is also evidence that some submergence-activated genes are clustered in the genome and that genes within these clusters are more highly expressed across angiosperms than submergence-activated genes outside of such clusters [30].
Alternative splicing (AS) offers a post-transcriptional mechanism for fine-tuning expression levels according to developmental stage, tissue or organ type, stress condition, or time [31]. In plants, the most common mode of AS is intron retention, whereas, in animals, exon skipping or shuffling is most frequent [32,33]. Most of these AS events (especially intron retention events) include premature stop codons, and transcripts are often subjected to degradation via nonsense-mediated decay, leading to downregulation of dosage [34][35][36]. Using RT-PCR, Syed et al. [37] found that in soybean, the SUB1 gene (a master controller of flooding resistance in rice) mediates flooding stress responses via AS, and PRR3 (a clock-associated gene) exhibits a flooding specific splicing pattern that appears to contribute to homeostasis. In maize, inbred lines tolerant of waterlogging exhibited differential splicing at prolyl 4-hydroxylases (which are thought to be oxygen sensors under hypoxia stress), compared to intolerant lines when exposed to flooding stress [38]. In Arabidopsis thaliana, transcriptomic analyses revealed water submergence-specific AS events in pathways regulating energy reserves such as gluconeogenesis [39]. However, we are unaware of other studies that have investigated genome-wide AS under flooding stress.
While transcriptomic responses to flooding stress have been investigated in a number of crops, much remains to be understood about which responses are common across many taxa and which are taxon-specific. In addition, previous studies typically have not examined interactions of genotypes with environment and time to determine the identity and timing of expression changes exclusive to tolerant lines. Here we report on transcriptomic responses to flooding in sunflower (Helianthus annuus L.), which is one of the world's most important oilseed crops, with global production valued at more than 20 Billion USD annually (FAOSTAT < http://www.fao.org/faostat/en/#data/QV>). However, yield is limited by several abiotic stresses, including drought, heat, salt, and flooding stress [40][41][42][43][44]. In North America, sunflowers sometimes experience short bouts of flooding in the spring, which results in the submergence of seeds, roots or even entire seedlings [45]. While such short-term flooding may not be lethal, yield and oil quality may be reduced [40].
Several potential adaptation strategies to flooding stress have been reported previously for both wild and cultivated sunflowers, including aerenchyma development and anaerobic energy production. Wample and Reid [46] reported that flooding induces the production of adventitious roots de novo, as well as swelling (hypertrophy) of the lower stem (hypocotyl) in cultivated sunflower. Both putative adaptations involve aerenchyma formation and are thought to represent a flooding escape strategy by promoting tissue gas exchange [46,47]. Likewise, evidence of increased anaerobic energy production has been reported for wild sunflowers that occupy habitats that frequently experience spring flooding [48]. A similar mechanism appears to be present in cultivated sunflower: flooding promotes the synthesis of ethanol, which is successfully metabolized thanks to an increase in alcohol dehydrogenase (ADH) activity [49].
Previously, we characterized phenotypic traits of 288 inbred sunflower lines from a cultivated sunflower association mapping (SAM) population under both flooding stress and benign conditions [44]. We measured multiple traits for possible flooding responses, such as biomass and chlorosis reduction, hypocotyl hypertrophy and adventitious root numbers, and identified candidate genes underlying these responses. We found that cultivated sunflower varied widely in flooding tolerance and that there was no inherent tradeoff between flooding tolerance and performance in the control treatment. On the other hand, we failed to find an effect of hypocotyl hypertrophy or adventitious root development on flooding tolerance, perhaps because they represent a general response to flooding across all cultivars. However, because of the ability of some of the genotypes to maintain growth without obvious negative impacts on photosynthesis (e.g., chlorosis), we suspected that increased anaerobic energy production might account for the success of tolerant phenotypes.
In the current study, we have extended our analyses of flooding tolerance in sunflower by studying expression and AS responses to flooding stress in one tolerant and susceptible line. We identified the genes that are differentially expressed or differentially spliced in response to flooding stress by generalized linear models or logistic regressions, respectively. Moreover, we used mixed model approaches to identify interactions between genotype, environment (in this case, control and flooding conditions), and time points during and post-flooding. We also conducted a gene co-expression network analysis to identify expression networks and to explore their connectedness and structure in resistant versus susceptible lines when responding to flooding stress. Our goal was not only to obtain a genome-wide assessment of molecular responses to flooding in sunflower but also to provide additional clues regarding how flooding tolerance is achieved. We hypothesize that pathways associated with aerenchyma formation and alcohol fermentation pathways play key roles in general and specific responses to flooding, respectively.

Plant Material and Flooding Stress Treatment
Out of 288 lines of the SAM population, we chose a resistant line (HA351) and a susceptible line (RHA428) based on vegetative growth performance and photosynthetic efficiency under the flooding stress ( Figure 1). The plants were germinated and grown following the protocol of Gao et al. [44]. Achenes were cleaned using 75% ethanol and 0.02% plant preservative mixture (PPM™; Plant Cell Technology Inc., Washington, DC, USA), scarified by cutting off a small segment from the blunt end of the achene, and then placed on damp filter paper in a Petri dish in the dark for three days. After germination, seed coats were removed, and seedlings were exposed to light for three days before they were transplanted into one-gallon pots. Each pot consisted of a 2:1 ratio of Sunshine Mix #1 (Sun Gro Horticulture, Vancouver, BC, Canada) and sand, respectively. The plants were grown in the greenhouse at the University of British Columbia (Vancouver, BC, Canada) under a regime of 14 h light and a temperature of 25 • C day/19 • C night. When the plants were three weeks old, half were placed in large plastic bins filled with tap water to impose flooding stress. The water level was kept at 1-2 cm above the soil surface throughout the stress treatment. After nine days, the pots were removed from the bins and returned to pre-flooding conditions for six days to allow recovery. Leaf and root tissue were collected from each plant, flash frozen in liquid nitrogen and kept at −80 • C until RNA extraction. Tissue was collected at three time points: (i) 24 h after flooding stress was imposed (1st day of flooding), (ii) the last day (8th day of flooding) of the flooding treatment and (iii) 6 days after the flooding treatment ended (6th day of recovery). These time points were chosen based on the results of Gao et al. [44].

RNA-Seq Data Processing
The quality of the 96 transcriptomes (75 bp paired-end raw reads, single-stranded) was assessed by FastQC (Version 0.11.8) (Andrews, Accessed 13 October 2017). Trimmomatic (Version 0.36) was used to discard low-quality reads and trim adaptor sequences [50]. The filtered reads were mapped to the sunflower reference genome [51] using STAR (Version 2.5.2b) with non-default parameters of alignIntronMax 10,000 [52]. We mapped reads against the reference genome instead of a transcriptome assembly to avoid multimapping of reads onto different isoforms of the same gene, which is critical for accurate analyses of alternative splicing. Mapped reads were quantified by featureCounts in the Subread package (Version 1.4.6) [53]. To further increase the effectiveness of the data analysis, HTSFilter (Version 1.12.0) [54] was used to remove genes with low constant expression levels. Analyses of root and leaf samples were conducted separately due to the strong difference between the two ( Figure S1).

Analyses of Differential Gene Expression and Alternative Splicing
For each time point, we identified differentially expressed genes (DEG) between the two conditions (control and flooding) using a false discovery rate (FDR) < 0.01 with fold- A total of 96 samples were collected: 2 genotypes × 2 tissues × 2 treatments × 3 time points × 4 biological replicates. For each biological replicate, leaflets and roots were pooled from three individual plants. Total RNA was extracted using the Ambion RNAqueous kit (AM1912, Invitrogen, Austin, TX, USA), and DNA contaminants were removed by the Ambion Turbo DNA-free kit (AM1907, Thermo Fisher, Waltham, MA, USA 02,451). A total of 96 libraries were prepared using the Illumina TruSeq stranded mRNA sample preparation kit and sequenced on an Illumina HiSeq 4000 (PE75) by Genome Quebec in Montréal (Québec, Canada). Sequence data from this project have been deposited in the National Center for Biotechnology Information Sequence Read Archive under accession number PRJNA492303.
A total of 48 seedlings (2 genotypes × 2 treatments × 3 time points × 4 biological replicates) were measured following the protocols of Gao et al. [44]: shoot biomass (SB, g of dry weight), leaf biomass (LB, g of dry weight), leaf area (LA, cm 2 ), leaf number (LN), plant height (PH, mm), stem diameter (SD, mm), chlorophyll concentration (CC, relative chlorophyll concentration unit via a chlorophyll meter), leaf mass per area (LMA, g/m 2 ), hypocotyl diameter (HD, mm) and adventitious roots development (number of adventitious roots, ARN) (Table S1). Raw trait values were then employed in a co-expression analysis (see below) to identify correlations with co-expression modules.

RNA-Seq Data Processing
The quality of the 96 transcriptomes (75 bp paired-end raw reads, single-stranded) was assessed by FastQC (Version 0.11.8) (Andrews, Accessed 13 October 2017). Trimmomatic (Version 0.36) was used to discard low-quality reads and trim adaptor sequences [50]. The filtered reads were mapped to the sunflower reference genome [51] using STAR (Version 2.5.2b) with non-default parameters of alignIntronMax 10,000 [52]. We mapped reads against the reference genome instead of a transcriptome assembly to avoid multi-mapping of reads onto different isoforms of the same gene, which is critical for accurate analyses of alternative splicing. Mapped reads were quantified by featureCounts in the Subread package (Version 1.4.6) [53]. To further increase the effectiveness of the data analysis, HTSFilter (Version 1.12.0) [54] was used to remove genes with low constant expression levels. Analyses of root and leaf samples were conducted separately due to the strong difference between the two ( Figure S1).

Analyses of Differential Gene Expression and Alternative Splicing
For each time point, we identified differentially expressed genes (DEG) between the two conditions (control and flooding) using a false discovery rate (FDR) < 0.01 with fold-change of −1 > log 2 FC or 1 < log 2 FC as calculated by edgeR (Version 3.18.1) [55], a Bioconductor package. Statistics of DEGs obtained from this approach are provided in Table S2. We converted raw read counts to reads per kilobase of transcript per million mapped reads (TPM) for downstream analyses. To visualize differences in expression, principal component analysis (PCA) was performed using prcomp in R package [56].
The alternative splicing analysis was conducted using custom Python scripts following the protocol of Tack et al. [57]. It identified the following major alternative splicing events: intron retention (IR), an alternative donor (ALTD), alternative acceptor (ALTA), alternative position (ALTP), exon skipping (SKIP), cryptic intron (CRIN) and cryptic exon (CREX). Reads mapping to the alternative and constitutive gene models were counted and quantified in terms of relative abundance per sample or percent splicing index (PSI) [58]. We employed logistic regression (FDR < 0.01 with fold-change of −1 > log 2 FC or 1 < log 2 FC) to identify differentially spliced (DAS) genes since the proportion of alternative splicing events is between 0 and 1 (calculated using a custom R script deposited on https://github.com/dejonggr/differential_as). We compared the log 2 FC between the two genotypes (pooling three time points together) to determine whether they differ in their response to flooding stress. Statistics for DAS genes obtained from this approach are provided in Table S3.

Mixed Effect Model Analyses of Variation in Gene Expression and Alternative Splicing
We used generalized linear mixed effect models to identify genes that responded differently in expression or AS to genotype, environment, time-point or the interaction between genotype and environment (G × E). We incorporated genotype, environment, time point and their interaction as fixed effects and individuals as random effects. We used a gamma family distribution with an inverse link function and a Gaussian family distribution with a logit transformation for expression (TPM value) and alternative splicing frequency data (PSI value), respectively. We used an FDR to correct for multiple comparisons across the entire experiment with an FDR cutoff of 0.1 to avoid false positives [59]. We performed independent contrasts to compare differences in expression or AS levels between genotypes, environment or time points and corrected for multiple comparisons using the false discovery rate. All genes were tested for G × E in the mixed-effects model.
In addition, we split each gene's data into three time points. Each time point was tested for G × E in a fixed-effects model. All analyses were done using the programming language R. Mixed effect models were performed using lme4 R package (Version 1.1-25) [60] and lmerTest (Version 3.1-3) [61] and analysis of variance was done using the "car" software package (Version 3.0-10) [62].

Gene Co-Expression Network Analysis
We used a weighted gene co-expression network analysis (WGCNA) (Version 1.69) [63] to identify co-expressed clusters of genes that respond to flooding stress in sunflower cultivars. The co-expression modules were assigned based on the TPM data of all expressed genes with a cutoff of at least three TPM in one condition throughout the four replicates (a total of 38,835 genes). Root and leaf samples were analyzed separately. Phenotypic trait data (Table S1) were used to test for correlations with the transcriptome data. The minimum module size was 100 genes, and the soft threshold power β was set to 12. The deep split was set to 2, and the network type was signed hybrid. Modules with the highest positive absolute correlation coefficient value with phenotypic data, and the lowest P value, were selected for further analyses. Genes with high connectivity were considered to be hub genes. Co-expression gene networks for the resistant and susceptible lines were generated independently using the same parameters to permit comparison. A consensus network analysis was conducted using the same parameters as described above to identify consensus modules and to better assess overall similarities and dissimilarities between the resistant and susceptible line. Networks were illustrated using Cytoscape (Version 3.7.2) [64].

Enrichment Comparisons
For all gene sets, gene ontology (GO) enrichment analyses were performed using agriGO Version 2.0 [65] with the Yekutieli multi-test adjustment method (false discovery rate < 0.1) [66], a minimum of five mapping entries and Plant GO slim annotations. The results are shown via REVIGO [67] TreeMaps that summarize enriched GO terms and their relationships.

RNA-Seq Data
To assess gene expression responses of cultivated sunflower to flooding stress, we conducted RNA-Seq on 96 libraries, representing two genotypes (resistant and susceptible), two environments (flooding and control), three time points (1st day of flooding treatment, 8th day of flooding treatment, and 6th day of recovery), two tissue types (root and shoot), and four biological replicates (Table 1). A total of 6,045,909,253 reads were obtained in this experiment (an average of 63 million reads per sample). After filtering out low-quality reads and genes with low constant expression levels, 3,016,817,944 clean Illumina RNA-Seq reads for 48,611 expressed genes were used to identify DEGs.

Differentially Expressed Genes (DEG) that Respond to Flooding Stress and Exhibit Genotype × Environment (G × E) Interactions
We found 2338 genes (1510 upregulated, 828 downregulated) from leaf tissue and 9089 genes (4417 upregulated, 4672 downregulated) from root tissue that were differentially expressed under flooding stress (E) (FDR < 0.01) ( Figure 2, Table S2). This root-dominant pattern of differential gene expression is shown in Table 2 and is displayed in Venn diagrams and scatter plots ( Figure S2). In addition, more differentially expressed genes were found during the flooding periods (time-point 1 and 2) compared to the recovery period (timepoint 3), especially in root tissue ( Figure 2 and Figure S2, Table 2). Interestingly, the number of differentially expressed genes in roots of the susceptible line at the time-point 3 (1247 genes) was far more than the resistant line (306 genes), suggesting that both genes responding to the flooding stress and those at the recovery stage, may contribute to flooding tolerance (Table 2 and Table S2, Figure S2). That is, DEGs from a mixed ANOVA of genotype (HA351 vs. RHA428) × treatment (control vs. flooding) interactions at time-point 1 in root tissue.

Differentially Expressed Genes (DEG) that Respond to Flooding Stress and Exhibit Genotype × Environment (G × E) Interactions
We found 2338 genes (1510 upregulated, 828 downregulated) from leaf tissue and 9089 genes (4417 upregulated, 4672 downregulated) from root tissue that were differentially expressed under flooding stress (E) (FDR < 0.01) ( Figure 2, Table S2). This root-dominant pattern of differential gene expression is shown in Table 2 and is displayed in Venn diagrams and scatter plots ( Figure S2). In addition, more differentially expressed genes were found during the flooding periods (time-point 1 and 2) compared to the recovery period (time-point 3), especially in root tissue (Figures 2 and S2, Table 2). Interestingly, the number of differentially expressed genes in roots of the susceptible line at the timepoint 3 (1247 genes) was far more than the resistant line (306 genes), suggesting that both genes responding to the flooding stress and those at the recovery stage, may contribute to flooding tolerance (Tables 2 and S2, Figure S2).   Principle components analysis (PCA) offered further insights into gene expression patterns ( Figure 3 and Figure S3). In root tissue, PC1 explained 23% of gene expression variation, which was partitioned mainly by treatment (control vs. flooding conditions). PC2 largely explained differences in genotype (resistant vs. susceptible lines), accounting for 13.4% of the variation, whereas PC3 appears to distinguish time-point 1 from the rest of the time-points and explains 10.4% of the variation. For leaf tissue, PC1 and PC2 explained 30.1% and 13.5% of gene expression variation, respectively, but the variation was not partitioned according to the fixed factors in our model. However, PC3 partitions samples by genotype (resistant versus susceptible) and explains 10.5% of the variation. This incongruence between leaf and root tissues may be accounted for by the direct contact of roots, but not leaves, to flooding stress. Indeed, the greater response of roots than leaves to flooding stress is seen throughout the experiment. Principle components analysis (PCA) offered further insights into gene expression patterns (Figures 3 and S3). In root tissue, PC1 explained 23% of gene expression variation, which was partitioned mainly by treatment (control vs. flooding conditions). PC2 largely explained differences in genotype (resistant vs. susceptible lines), accounting for 13.4% of the variation, whereas PC3 appears to distinguish time-point 1 from the rest of the timepoints and explains 10.4% of the variation. For leaf tissue, PC1 and PC2 explained 30.1% and 13.5% of gene expression variation, respectively, but the variation was not partitioned according to the fixed factors in our model. However, PC3 partitions samples by genotype (resistant versus susceptible) and explains 10.5% of the variation. This incongruence between leaf and root tissues may be accounted for by the direct contact of roots, but not leaves, to flooding stress. Indeed, the greater response of roots than leaves to flooding stress is seen throughout the experiment. To investigate which genes show differential expression between the resistant and susceptible lines due to the flooding treatment and/or the timing of the treatment, we used To investigate which genes show differential expression between the resistant and susceptible lines due to the flooding treatment and/or the timing of the treatment, we used a mixed model approach ( Figure 4; Table S2). After the FDR correction (0.1 or less), we found that G × E interactions were significant for 1256 genes in leaf tissue and 1093 genes in root tissue. Of these genes, we then searched for the subset that differed significantly (twofold, FDR < 0.1) between lines under flooding stress but not under control conditions. In leaf tissue (G × E_leaf), we found 8 genes that were significantly upregulated and 31 genes that were significantly downregulated in the resistant line compared to the susceptible line. In the root tissue (G × E root), we found 20 genes that were significantly upregulated and 41 genes that were significantly downregulated in the resistant line relative to the susceptible line. These genes represent good candidates for explaining the differences between the two lines in their resistance to the flooding stress (Table S2-G × E leaf and G × E root).
genes that were significantly downregulated in the resistant line compared to the susceptible line. In the root tissue (G × E root), we found 20 genes that were significantly upregulated and 41 genes that were significantly downregulated in the resistant line relative to the susceptible line. These genes represent good candidates for explaining the differences between the two lines in their resistance to the flooding stress (Table S2-G × E leaf and G × E root). For the time-point analysis, we searched for the genes that were significantly different between lines at one or more time-points in the flooding treatment but not in the control treatment (Table S2). In leaf tissue, we found 91 genes that were upregulated and 178 genes that were downregulated in the resistant versus susceptible line at time-point 1 (G × E1_leaf). At time-points 2 and 3 (G × E2-3_leaf), 82 and 59 genes were upregulated, and 182 and 208 genes were downregulated, respectively, in the resistant versus susceptible line. In root tissue, we found 121, 1,101, and 170 genes that were upregulated and 282, 287, and 118 genes that were downregulated at time-points 1, 2, and 3, respectively, in the resistant versus susceptible line (G × E1-3_root). The identity and function of these genes may offer clues regarding how shifts in the timing of expression of particular categories of genes contribute to flooding resistance (Table S2-G × E1-3_leaf and G × E1-3_root).

Differentially Spliced Genes (DAS) that Respond to the Flooding Stress
We identified AS events and estimated their frequencies using the transcriptome data described above (Table S4). As expected, intron retention events were most frequent and exon skipping least frequent in both tissues. All AS events included premature stop codons, suggesting that they could be subject to nonsense-mediated mRNA decay (NMD) and downregulate the overall expression of the gene. An AS ratio was calculated to determine the frequency of reads that mapped to the alternatively spliced gene model compared to the constitutive gene model. Thus, upregulation or downregulation of AS refers to an increase or decrease, respectively, of the AS ratio. We found 892 AS events (132 upregulated, 760 downregulated) from leaf tissue and 2734 events (1578 upregulated, 1156  For the time-point analysis, we searched for the genes that were significantly different between lines at one or more time-points in the flooding treatment but not in the control treatment (Table S2). In leaf tissue, we found 91 genes that were upregulated and 178 genes that were downregulated in the resistant versus susceptible line at time-point 1 (G × E1_leaf). At time-points 2 and 3 (G × E2-3_leaf), 82 and 59 genes were upregulated, and 182 and 208 genes were downregulated, respectively, in the resistant versus susceptible line. In root tissue, we found 121, 1101, and 170 genes that were upregulated and 282, 287, and 118 genes that were downregulated at time-points 1, 2, and 3, respectively, in the resistant versus susceptible line (G × E1-3_root). The identity and function of these genes may offer clues regarding how shifts in the timing of expression of particular categories of genes contribute to flooding resistance (Table S2-G × E1-3_leaf and G × E1-3_root).

Differentially Spliced Genes (DAS) That Respond to the Flooding Stress
We identified AS events and estimated their frequencies using the transcriptome data described above (Table S4). As expected, intron retention events were most frequent and exon skipping least frequent in both tissues. All AS events included premature stop codons, suggesting that they could be subject to nonsense-mediated mRNA decay (NMD) and downregulate the overall expression of the gene. An AS ratio was calculated to determine the frequency of reads that mapped to the alternatively spliced gene model compared to the constitutive gene model. Thus, upregulation or downregulation of AS refers to an increase or decrease, respectively, of the AS ratio. We found 892 AS events (132 upregulated, 760 downregulated) from leaf tissue and 2734 events (1578 upregulated, 1156 downregulated) from root tissue that varied significantly in frequency under flooding stress (Table S3). Similar to the DEGs described above, differentially spliced events were more than twice as common in root compared to leaf tissue. A total of 69 genes (11 upregulated, 58 downregulated) from leaf tissue and 97 genes (56 upregulated, 41 downregulated) from root tissue were differentially spliced in both genotypes under flooding stress.
We also conducted a PCA of differentially spliced genes under flooding stress ( Figure S4). Unlike the expression data, PC1 did not differentiate DAS in any variable. On the other hand, PC2 of leaf tissue (7.6% of variation) and PC3 of root tissue (8.4% of variation) were associated with the genotypic differences.

Putative Functions of Differentially Expressed Genes
Of the genes that were differentially expressed under flooding stress (E leaf), those upregulated in leaf tissue were enriched in 71 GO terms (Table S5). Many of the significant GO terms were related to protein translation and response to an extracellular stimulus (biological process; Figure S5A), ribosome (cellular component, Figure S5B), and structural molecule activity (molecular function; Figure S5C). Fifty-five GO terms were significantly enriched for genes that were downregulated in leaf tissue (E_leaf); GO terms relating to secondary metabolism, homeostatic process and photosynthesis (biological process; Figure S5D), cell wall and membrane (cellular component, Figure S5E), and catalytic activity and binding (molecular function; Figure S5F) were especially frequent. Overall, these results indicate that in leaf tissue under flooding stress, genes involved in the production of extracellular stimulus-responsive peptides (possibly due to ethylene and ROS) were upregulated, whereas those related to secondary metabolism, energy production or photosynthesis were downregulated. Slowing down these and other metabolic pathways may contribute to flooding tolerance.
For root tissue under flooding stress, upregulated genes (E_root) were enriched in 75 GO terms (Table S5). The most significant GO terms were enriched in response to an extracellular stimulus (biological process; Figure S5G), cell wall and membrane (cellular component, Figure S5H), catalytic activity and transferase activity (molecular function; Figure S5I). Genes that were significantly downregulated in root tissue (E_root) were enriched in 84 GO terms, especially secondary metabolism and response to abiotic stress (biological process; Figure S5J), cell wall and membrane (cellular component, Figure S5K), catalytic activity and transferase activity (molecular function; Figure S5L). Overall, root tissue responded to flooding stress as expected, by upregulating expression levels of flooding responsive genes. In addition, many metabolism-related activities appear to have slowed, presumably due to a lack of oxygen and energy.
We also conducted GO enrichment analyses for the DEGs identified in the mixed models (Table S5) In the leaf tissue (G × E_leaf), only eight genes with significant G × E were upregulated, which is too few for GO enrichment analysis. Nonetheless, several of the genes were notable and include HanXRQChr08g0214611 (zinc finger, RING/FYVE/PHD-type known as XERICO) and HanXRQChr14g0450491 (high-affinity nitrate transporter 2.7). Thirty-one genes with significant G × E interactions were downregulated; this set of genes was enriched in response to stress (biological process), extracellular region (cellular component), and catalytic activity and binding (molecular function). Results of G × E analyses were similar for each time point (G × E1-3_leaf): signaling-related genes were upregulated, and response to stress genes were downregulated. This implies that stress response genes in the resistant line are less responsive to flooding than in the susceptible line, which may underlie the greater growth performance of the former.
In the root tissue (G × E_root), upregulated genes exhibiting significant G × E interactions were related to GO: 0003824 "catalytic activity" (FDR = 0.031; molecular function) and GO: 0005576 "extracellular region" (FDR = 0.049; cellular component), whereas, downregulated genes were enriched in cellular component GO term of "intracellular organelle part" (FDR = 0.0025). Genes with significant G × E interactions at time-point 1 (G × E1_root) were most strongly enriched in GO:0003824 "catalytic activity" for both up-and downregulated genes; significant genes at time-point 2 (G × E2_root) were enriched in GO: 0009628 "response to abiotic stimulus" (FDR = 0.085) in the upregulated fraction and GO: 0003824 "catalytic activity" (FDR = 8.00E-08) in the downregulated fraction, and significant genes at time-point 3 (G × E3_root) were enriched in GO:0050896 "response to stimulus" (FDR = 3.60E-05) in the upregulated subset and GO:0003824 "catalytic activity" (FDR = 0.001) in the downregulated subset. These enrichment patterns indicate that the resistant and susceptible line differ significantly both in the kinds of genes that respond to flooding stress and the timing of the response. Genes involved in catalytic activity appear to play an especially critical role in flooding response. Notably, four ADH genes (HanXRQChr05g0161031, HanXRQChr05g0161051, HanXRQChr06g0179461 and HanXRQChr09g0243891) were upregulated at time-point 1 (G × E1_root), as well as two additional ADH genes (HanXRQChr05g0161051 and HanXRQChr17g0559391) at timepoint 2 (G × E2_root). This implies that rapid upregulation of the fermentation pathway contributes to the greater flooding tolerance of the resistant line.

Putative Function of Differentially Spliced Genes
Among the genes that had a significant change in AS frequency under flooding (E_leaf), those with increased AS in leaf tissue were significantly enriched in carbohydrate metabolism (biological process; Figure S6A), intracellular organelle part (cellular component, Figure S6B), and pyrophosphatase activity (molecular function; Figure S6C). Genes with decreased AS in leaf tissue were significantly enriched in embryo development (biological process; Figure S6D), intracellular organelle part (cellular component, Figure S6E), and in transferase activity (molecular function; Figure S6F). Increased AS frequency may reduce the expression level of metabolic genes via NMD, further downregulating metabolic pathways under flooding stress.
Genes with significant increases in AS frequencies under flooding stress in root tissue (E_root) were enriched in translation (biological process; Figure S6G), cell part (cellular component, Figure S6H), and nuclease activity (molecular function; Figure S6I). Genes with decreases in AS were enriched in embryo development (biological process; Figure S6J), cell part (cellular component, Figure S6K), and nuclease activity (molecular function; Figure S6L). This pattern of increased AS of metabolic process and decreased AS of embryo development are accord with those in leaf tissue, suggesting the regulation of AS in both leaf and root tissue work in a similar way and likely acts to fine-tune gene expression responses to environmental change.
The mixed model analysis (Table S3) indicated that very few AS events exhibit significant G × E interactions, so we were unable to be able to conduct GO enrichment analyses. Therefore, individual DAS events were explored. For example, AS frequencies of HanXRQChr09g0273081 (elongation factor 1-delta) and HanXRQChr11g0325531 (translation initiation factor IF2/IF5) were reduced at time-point 3 in leaf tissue (G × E3_leaf), consistent with the reactivation of cell expansion and growth in the recovery stage. In root tissue, HanXRQChr05g0161051 (alcohol dehydrogenase) had reduced AS in the resistant line at time-point 1 (G × E1_root), and this may contribute to more efficient expression via constitutive splicing in the fermentation pathway. HanXRQChr12g0371921 (myctype, basic helix-loop-helix) showed an increase in AS in the resistant line at time-point 1 (G × E1_root), and its molecular function (GO: 0003700) is related to DNA-binding transcription factor activity.
Overall, the changes in AS reinforce expression changes under flooding stress, with strong negative correlations observed (r ranges from −0.57 to −0.61 and P from 5.01 × 10 −11 to 1.41 × 10 −49 ) regardless of tissue type or treatment ( Figure S11).

Gene Co-Expression Network Modules Underlying Responses to Flooding Stress
Co-expression networks were constructed using all 96 RNA-Seq transcriptomes, but root and leaf transcripts were analyzed separately due to the large differences between the two tissue types. We generated network modules for resistant and susceptible lines independently using standard parameters (see Materials and Methods) and compared them via consensus analysis following the protocol of Langfelder and Horvath [63]. We were particularly interested in correlations that differed significantly between resistant and susceptible data sets in strength or sign, thereby permitting identification of unique gene co-expression networks in the flooding resistant sunflower cultivar.
The leaf tissue network for the resistant line was comprised of 34 clusters of coexpressed genes. These module eigengenes (ME) were identified by the first component in a PCA, which was generated by the WGCNA package. Module-trait associations are presented by a heat map ( Figure S7A), which illustrates the correlations between the MEs and phenotypic traits. Chlorophyll concentration is closely connected to flooding resistance in sunflower due to its impact on photosynthesis [44]. We observed that genes in module ME-black (Table S7) in the heat map are most positively correlated (correlation coefficient = 0.57, P = 0.004) with chlorophyll concentration. ME-black consists of genes that were enriched in the following GO terms (Table S8): response to an extracellular stimulus (biological process; Figure S7B), plasma membrane (cellular component, Figure S7C), and kinase activity (molecular function; Figure S7D). The ME-salmon module from the consensus data analysis, which shows substantial overlap with the ME-black module from the resistant line ( Figure S8A), has a significant positive correlation (correlation coefficient = 0.53, P = 0.007) with chlorophyll concentration in the resistant line ( Figure S8B), but a non-significant negative correlation in the susceptible line (correlation coefficient = −0.28, P = 0.2; Figure S8C). Thus, these gene co-expression modules ( Figure S8D) may contribute to the different phenotypic responses to flooding between the resistant and susceptible line.
We further explored the ME-black module to determine whether it includes the differentially expressed genes identified by the mixed model, and if so, how they are connected or embedded within the gene co-expression network. We visualized the network by Cytoscape, only including the top 10 hub genes (with the highest degree of connection) and the differentially expressed genes in the G × E_leaf analysis (Table S7; Figure S7E). Many of the hub genes were transcription factors, including HanXRQChr01g0029281 (zinc finger, PHD-type), HanXRQChr04g0103501 (WRKY domain; Zn-cluster domain), HanXRQChr05g0160581 (NAC domain transcriptional regulator), and HanXRQChr07g0191711 (zinc finger, RING/ FYVE/PHD-type). Fifteen differentially expressed genes (nine upregulated and six downregulated) identified in the G × E analysis at time-point 1 (G × E1_leaf) were also found (Table S2) in the ME-black module. Interestingly, HanXRQChr13g0404381 (lipoxygenase; lipoxygenase, domain 3) was upregulated, and HanXRQChr09g0252851 (lipase class 3related protein) were downregulated in the resistant compared to the susceptible line, suggesting regulation of lipids is critical at the beginning of flooding stress. The latter gene continued to be downregulated in the resistant line at time-point 2 (G × E2_leaf) (Table S2). By time point three (G × E3_leaf), all 15 differentially expressed genes were downregulated in the resistant versus susceptible line (Table S2), consistent with our hypothesis that the resistant line recovers more quickly than the susceptible line by turning down stress response pathways.
In the root tissue, we identified 38 clusters of co-expressed genes in the resistant line. The module-trait associations are presented by a heat map ( Figure S9A). Since we do not have root phenotypic trait data, any correlations between root tissue expression and phenotypic changes are expected to be indirect (e.g., development or growth rate). We chose to focus on the ME-red module (Table S7) since it has the strongest positive correlation among module-trait relationships, and it is positively correlated with leaf biomass (correlation coefficient = 0.9, P = 3 × 10 −9 ) and leaf area (correlation coefficient = 0.72, P = 8 × 10 −5 ), both which are significantly reduced by flooding stress (Table S1; also see [44]). ME-red consists of the genes that were enriched in GO terms (Table S8) related to the response to endogenous stimulus and carbohydrate metabolism (biological process; Figure S9B), cell wall (cellular component, Figure S9C), and catalytic activity (molecular function; Figure S9D). The ME-red module overlapped with a resistant-susceptible consensus module (ME-purple) ( Figure S10A). ME-purple of the consensus data have a significant correlation (correlation coefficient = 0.55, P = 0.005) with the leaf area in the resistant line ( Figure S10B), but not in the susceptible line (correlation coefficient = −0.04, P = 0.9; Figure S1C). The lack of consensus for this trait between resistant and susceptible lines is indicated as NA in Figure S10D.
The ME-red module in the root tissue of the resistant line was further analyzed by identification of hub genes. As before, we used Cytoscape for network visualization with the top 10 hub genes (with the highest degree of connection), as well as differentially expressed genes from the G × E analysis (Table S7; Figure 5). The hub genes with the highest degree of connectivity include HanXRQChr09g0248141 (hypoxia-responsive family protein), HanXRQChr12g0380481 (transcription factor TGA like domain), and HanXRQChr15g0464051 (Thiamine pyrophosphate dependent pyruvate decarboxylase family protein). Twenty-one differentially expressed genes exhibiting significant G × E interactions at time-point 1 (G × E1_root) were all upregulated in the resistant line (Table S2) and belong to the ME-red module. Among the DEGs, three ADH homologs (HanXRQChr05g0161031, HanXRQChr05g0161051, HanXRQChr06g0179461) were found. One of these (HanXRQChr05g0161051) was also upregulated at time-point 2 (G × E2_root) (Table S2). This is consistent with the results from the mixed model analyses described previously, suggesting that more rapid and stronger upregulation of the fermentation pathway in the resistant line may contribute to its greater tolerance to flooding. S10B), but not in the susceptible line (correlation coefficient = −0.04, P = 0.9; Figure S1C). The lack of consensus for this trait between resistant and susceptible lines is indicated as NA in Figure S10D.
The ME-red module in the root tissue of the resistant line was further analyzed by identification of hub genes. As before, we used Cytoscape for network visualization with the top 10 hub genes (with the highest degree of connection), as well as differentially expressed genes from the G × E analysis (Table S7; Figure 5). The hub genes with the highest degree of connectivity include HanXRQChr09g0248141 (hypoxia-responsive family protein), HanXRQChr12g0380481 (transcription factor TGA like domain), and HanXRQChr15g0464051 (Thiamine pyrophosphate dependent pyruvate decarboxylase family protein). Twenty-one differentially expressed genes exhibiting significant G × E interactions at time-point 1 (G × E1_root) were all upregulated in the resistant line (Table  S2) and belong to the ME-red module. Among the DEGs, three ADH homologs (HanXRQChr05g0161031, HanXRQChr05g0161051, HanXRQChr06g0179461) were found. One of these (HanXRQChr05g0161051) was also upregulated at time-point 2 (G × E2_root) (Table S2). This is consistent with the results from the mixed model analyses described previously, suggesting that more rapid and stronger upregulation of the fermentation pathway in the resistant line may contribute to its greater tolerance to flooding.  Venn diagram of ME-red and differentially expressed genes with significant two-way genotype by treatment (i.e., environment) effects. GE indicates genotype × environment interactions, and the number indicates the time-point. (B) Visualization of network connections among the most connected genes in ME-red, as generated by Cytoscape software. Orange circles are hub genes with the highest connectivity. Green circles are alcohol dehydrogenase (ADH) homologs. White circles are differentially expressed genes with significant genotype by treatment interactions when combined across the three time points (except for the ADH homologs).

Comparison with Evolutionarily Conserved Submergence Activated Genes
A comprehensive study of submergence upregulated gene families (SURFs) was reported recently by Reynoso et al. [30]. They analyzed root tissue of four flowering plant species: rice, legume (Medicago truncatula), domesticated tomato (Solanum lycopersicum) and its dryland-adapted wild relative (Solanum pennellii) to understand variation and conservation of SURFs among species.
In our root expression data, we identified 35 homologs of conserved SURFs that responded significantly to flooding, of which 28 were upregulated and seven downregulated ( Table S2). Two of these SURFs (AT1G77120-alcohol dehydrogenase and AT4G33070pyruvate decarboxylase 1) were found to upregulated in the resistant line relative to the susceptible line at time-points 1 and 2, and at time-point 1, respectively (Table S2). Pyruvate decarboxylase initiates the ethanol fermentation process under anaerobic conditions by converting pyruvate into acetaldehyde and carbon dioxide, providing additional evidence for the importance of accelerated upregulation of the alcohol fermentation pathway.
With respect to alternative splicing, a significant increase in AS frequency was found for six homologs of SURF genes, whereas for three genes, a significant decrease in AS frequency was observed (Table S3). Notably, one of the genes with reduced AS was HanXRQChr05g0161021 (AT1G77120-alcohol dehydrogenase), further increasing the number of functional transcripts of this gene under flooding stress.

Case Study of the Alcohol Fermentation Pathway
Given previous evidence pointing towards the importance of the alcohol fermentation pathway in flooding tolerance in sunflower [49], as well as differences between the resistant and susceptible line observed in the present study, we decided to take a more comprehensive look at this gene family. We retrieved 65 ADH homologs in the sunflower genome (Table S9) [51]. In leaf tissue, nine ADH homologs were differentially expressed (mostly upregulated) under flooding stress; however, none had significant G × E interactions ( Table S9), suggesting that upregulation of ADH homologs is a general response to flooding stress in leaf tissue. In contrast, many ADH homologs were differentially expressed under flooding stress (both up and down) in the root tissue (Table S9), including five of which were upregulated in the resistant line, but not the susceptible line. Four of these were upregulated at time-point 1 (G × E1_root) (HanXRQChr05g0161031, HanXRQChr05g0161051, HanXRQChr06g0179461, and HanXRQChr09g0243891) (Table S9), and three were included in the red module of the resistant line that was described in the co-expression network analysis.
Noticeably, four ADH homologs that were responsive to flooding appear to be tandemly duplicated on chromosome 5 ( Figure 6A). All four had higher expression in the resistant versus susceptible line under the flooding stress, with the strongest difference observed for HanXRQChr05g0161051 (G × E1_root; FDR = 6.87 × 10 −12 ). In addition, as mentioned above, this gene showed a reduction in AS frequency in the resistant line relative to the sensitive line at time-point 1 (G × E1_root; FDR = 9.95 × 10 −6 ; Figure 6B).

Root Dominant Responses to Flooding Stress in Sunflower
Changes in expression and AS were greater in root than shoot tissue in all analyses and at all time-points. This root-dominant response to flooding has been observed in several other taxa, including Taxodium, a flooding-tolerant conifer. Taxodium plantlets exhibited root dominant differential expression of the antioxidative defense, glycolysis, and fermentation pathways, which play a key role in adaptive flooding responses by reducing damage from ROS species in the former pathway and maintaining ATP production in the latter two pathways [68]. Likewise, gray poplar (Populus × canescens) displayed greater transcriptional change in roots than leaves under flooding stress and upregulated gene pathways were similar to those reported for Taxodium [69]. However, there are exceptions to this pattern. For example, van Veen et al. [39] reported greater expression divergence in shoot than root tissue in Arabidopsis under flooding stress. This may be a consequence

Root Dominant Responses to Flooding Stress in Sunflower
Changes in expression and AS were greater in root than shoot tissue in all analyses and at all time-points. This root-dominant response to flooding has been observed in several other taxa, including Taxodium, a flooding-tolerant conifer. Taxodium plantlets exhibited root dominant differential expression of the antioxidative defense, glycolysis, and fermentation pathways, which play a key role in adaptive flooding responses by reducing damage from ROS species in the former pathway and maintaining ATP production in the latter two pathways [68]. Likewise, gray poplar (Populus × canescens) displayed greater transcriptional change in roots than leaves under flooding stress and upregulated gene pathways were similar to those reported for Taxodium [69]. However, there are exceptions to this pattern. For example, van Veen et al. [39] reported greater expression divergence in shoot than root tissue in Arabidopsis under flooding stress. This may be a consequence of different resistance strategies, with Arabidopsis exhibiting rapid shoot growth to escape submergence, whereas sunflower (for example) appears to upregulate anaerobic energy production in roots and produces adventitious roots. In contrast to shootdominant expression divergence, Arabidopsis shows greater differential AS in roots when exposed to flooding stress [39], an observation that aligns with our study. Establishing root dominant AS as a general rule, however, awaits other studies of AS responses to flooding across multiple tissue types.

General Responses to Flooding
We observed two main responses to flooding stress (E_leaf and root) that were common to both the resistant and susceptible lines: upregulation of ethylene-responsive transcription factors and downregulation of metabolic, energy production, and photosynthesisrelated pathways. The accumulation of ethylene gas is considered to be the most consistent, timely and reliable indicator of early flooding stress [70], and ethylene-responsive transcription factors (ERFs) play an early and critical role in mediating transcriptomic responses to flooding stress [71]. Ethylene regulates a hormonal cascade (abscisic acid, gibberellic acid and auxin) that along with ROS promote adventitious root growth [11], and several previous studies have shown that upregulation of ERFs is positively correlated with adventitious root formation (e.g., in cucumber [72] and rice [73]). ERFs were significantly upregulated under flooding stress in sunflower at the first time point (E_leaf and root) and may be responsible for the increase in adventitious rooting seen under flooding stress in our phenotypic data (Table S1), as well as in our previous large-scale GWA study of sunflower cultivars [44]. It is important to keep in mind, though, that upregulation of ERFs and increased adventitious rooting are seen in both the resistant and susceptible lines, so this appears to be a general response to flooding and not specific to the resistant genotype.
Another common response to flooding is reduced photosynthesis activity due to lower availability of oxygen, which can lead to chlorosis, leaf abscission and a decline in metabolic rates [12,20,74,75]. Better control of metabolic pathways and energy use can reduce negative impacts from flooding stress [15,76]. In sunflower, we observed phenotypic evidence of reduced photosynthesis activity ( Figure 1, Table S1), and similar results were reported previously by Gao et al. [44]. Reduced photosynthetic activity can also be inferred from the downregulation of metabolic, energy production, and photosynthesisrelated pathways in both the resistant and susceptible cultivars in the flooding treatment ( Figure S5D). The downregulation of photosynthesis-related gene pathways following flooding has previously been reported in several other plant species, including soybean [27] and canola [77].

Changes in AS Levels Reinforce Expression Differences
Alternative splicing is a post-transcriptional mechanism that can increase the diversity of the proteome and regulate gene expression [78][79][80]. The latter typically occurs via the expression of alternative transcripts that are targeted to NMD [81,82]. Increases in the frequency of AS for a given gene effectively lower the expression of functional transcripts. This was observed by Filichkin and Mockler [83], who reported an increase in nonsense isoforms (and a lower proportion of functional transcripts) in Arabidopsis thaliana following temperature stress. On the other hand, decreasing frequencies of AS can restore the function of genes that are inhibited by AS. For example, the DEHYDRATION-RESPONSIVE ELEMENT BINDING 2 (DREB2) transcription factor produces a non-functional transcript under normal conditions but a functional transcript in response to abiotic stress [84,85].
In our experiment, regardless of genotype, tissue type or treatment, levels of AS were strongly negatively correlated with levels of gene expression ( Figure S11), thereby reinforcing expression divergence. For example, metabolism-related gene pathways exhibited decreased expression levels and increased AS frequencies (i.e., fewer functional transcripts) under flooding stress (E_leaf) (Figures S5 and S6). While the number of genes that had significant differences in both expression and AS was relatively small, in most cases, differences were in the opposite direction, leading to reinforcing effects on the number of functional transcripts. This was most pronounced in leaf tissue (E_leaf), where genes that showed increased expression and downregulated AS were most frequent ( Figure S12). Such a complementary pattern between expression and AS levels has been previously reported in a study of homologous expression and AS patterns in canola under several abiotic stresses [86]. A focus of future studies should be to validate these inferences using proteomic approaches.

Expression Changes Specific to the Flooding Resistant Cultivar
The flooding resistant cultivar differed from the susceptible cultivar in three main ways: earlier and stronger upregulation of the alcohol fermentation pathway muted expression of stress-responsive genes in leaf tissue and more rapid recovery of pre-stress expression levels. Below we discuss how these differences may lead to greater resistance or tolerance.
The alcohol fermentation pathway represents a critical metabolic, and physiological response mechanism to hypoxia in higher plants by providing the energy (NAD + ) needed for glycolysis to continue in the absence of oxygen [7][8][9]87]. The pathway is a highly conserved adaptation to low oxygen, and upregulation of the pathway can enhance a plant's tolerance to flooding stress [10,88].
Alcohol dehydrogenase (ADH), the core enzyme in the alcohol fermentation pathway, is considered to play an especially important role in flooding tolerance because it both detoxifies toxic acetaldehyde and regenerates oxidized NAD + [87]. We explored the expression of ADH homologs (orthologous to AT1G77120 in Arabidopsis thaliana). Our mixed model analysis of genotype x treatment across the three different time points revealed that a number of ADH homologs showed earlier and/or stronger expression in the resistant than susceptible genotype when under flooding stress (G × E1_root vs. G × E2_root or G × E3_root). These differences were further reinforced by significant downregulation of ADH AS levels (leading to a larger fraction of functional transcripts) in the resistant line relative to the susceptible line. Notably, five ADH homologs on chromosome 5 appear to be the product of tandem duplication, and they were upregulated much more strongly under flooding in the resistant than susceptible line (G × E1_root) ( Figure 6A). Such duplications offer a means for increasing expression levels [89] and are known to contribute to abiotic stress tolerance in plants. A well-known example is the polygenic SUBMERGENCE 1 (SUB1) locus, which confers tolerance to flooding in rice [90]. In addition, most ADH homologs were found in a flooding resistance-related co-expression network in root tissue, further supporting the importance of this gene family. The genes found in this co-expression network largely overlap with hypoxia resistance networks previously reported for other plants [39,91], except for the inclusion of alcohol fermentation pathway in the sunflower network. The hub genes in this co-expression network, such as HanXRQChr08g0216771 (phosphofructokinase) and HanXRQChr12g0380481 (transcription factor TGA), likely represent key regulators of the alcohol fermentation pathway, although this claim awaits functional validation.
In addition to the earlier and stronger expression of the alcohol fermentation pathway, our results indicate that the upregulation of stress-responsive genes is muted in leaf tissue of the resistant genotype and that the resistant genotype recovers pre-stress expression levels more quickly than the sensitive genotype. The rapid recovery of metabolic homeostasis is often associated with stress tolerance [92,93], but whether this is a cause or consequence of increased resistance is not clear in most cases. In sunflower, it may be that the sensitive genotype over-responds to flooding, leading to reduced growth rates. This has been reported for rice, in which overexpression of a stress-responsive gene, OsHOX24, reduces root and shoot growth [94]. Alternatively, the greater energy production in roots of the resistant genotype (and/or other changes not discussed here) may reduce the negative impacts of flooding, leading to a less severe response and more rapid recovery.

Conclusions
We conducted a large-scale transcriptomics study to investigate the eco-physiological responses of sunflower, an important oilseed crop, to flooding stress. Using mixed model and co-expression network approaches, we explored the differential expression and AS patterns to identify genes and genetic pathways that underlie general responses to flooding stress, as well as those specific to a resistant cultivar. We found that upregulation of ethylene-responsive transcription factors and adventitious root growth occurred in both the sensitive and resistant cultivated lines and appear to represent a general critical response to flooding stress. Our hypothesis that differential expression and AS of the alcohol fermentation pathway contributes to the greater tolerance of the resistant line was supported. Greater flooding tolerance of the resistant line is also associated with muted expression of stress-responsive genes in leaf tissue and with a more rapid recovery of "normal" expression after the stress period is over. A goal for the future will be to extend our results across more genotypes and to validate expression and AS patterns with proteomic data.  Figure S7: Gene co-expression network analysis of leaf tissue of resistant line; Figure S8: Consensus modules in leaf tissue; Figure S9: Gene co-expression network analysis of root tissue of resistant line; Figure S10: Consensus modules in root tissue; Figure S11: Scatter plot of Pearson's product moment correlation coefficient between differential expression (DEG) and differential splicing (DAS); Figure S12: Overlap between differentially expressed genes (DEGs) and differentially spliced genes (DAS) under flooding, Table S1: Phenotypic measurement data in this study; Table S2: Differentially expressed gene responding to different variables; Table S3: Differentially spliced gene responding to different variables; Table S4: Frequency of alternative splicing events estimated from the sunflower transcriptome data reported here; Table S5: Gene ontology enrichment analysis results of differentially expressed genes; Table S6: Gene ontology enrichment analysis results of differentially spliced genes; Table S7: Gene co-expression network modules; Table S8: Gene ontology enrichment analysis results of gene co-expression network modules; Table S9. Alcohol dehydrogenase homologs in sunflower.