Science, Faculty of
Non UBC
Botany, Department of
DSpace
Agronomy 11 (1): 92 (2021)
University of British Columbia. Biodiversity Research Centre
Lee, Joon Seon
Gao, Lexuan
Guzman, Laura Melissa
Rieseberg, Loren H.
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.
agronomyArticleGenome-Wide Expression and Alternative Splicing inDomesticated Sunflowers (Helianthus annuus L.) underFlooding StressJoon Seon Lee 1,* , Lexuan Gao 2 , Laura Melissa Guzman 3,4 and Loren H. Rieseberg 1Citation: Lee, J.S.; Gao, L.; Guzman,L.M.; Rieseberg, L.H. Genome-WideExpression and Alternative Splicingin Domesticated Sunflowers(Helianthus annuus L.) under FloodingStress. Agronomy 2021, 11, 92.https://doi.org/10.3390/agronomy11010092Received: 4 November 2020Accepted: 2 January 2021Published: 6 January 2021Publisher’s Note: MDPI stays neu-tral with regard to jurisdictional clai-ms in published maps and institutio-nal affiliations.Copyright: © 2021 by the authors. Li-censee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms and con-ditions of the Creative Commons At-tribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).1 Department of Botany and Biodiversity Research Centre, University of British Columbia,Vancouver, BC V6T 1Z4, Canada; lriesebe@mail.ubc.ca2 CAS Center for Excellence in Molecular Plant Sciences, Institute of Plant Physiology and Ecology,Chinese Academy of Sciences, Shanghai 200032, China; lexuan@sippe.ac.cn3 Evolutionary and Behavioural Ecology Research Group, Department of Biological Sciences,Simon Fraser University, Burnaby, BC V5A 1S6, Canada; guzman@zoology.ubc.ca4 Department of Zoology and Biodiversity Research Centre, University of British Columbia,Vancouver, BC V6T 1Z4, Canada* Correspondence: joon.lee@botany.ubc.ca; Tel.: +1-604-827-3535Abstract: Approximately 10% of agricultural land is subject to periodic flooding, which reducesthe growth, survivorship, and yield of most crops, reinforcing the need to understand and enhanceflooding resistance in our crops. Here, we generated RNA-Seq data from leaf and root tissue ofdomesticated sunflower to explore differences in gene expression and alternative splicing (AS)between a resistant and susceptible cultivar under both flooding and control conditions and atthree time points. Using a combination of mixed model and gene co-expression analyses, we wereable to separate general responses of sunflower to flooding stress from those that contribute to thegreater tolerance of the resistant line. Both cultivars responded to flooding stress by upregulatingexpression levels of known submergence responsive genes, such as alcohol dehydrogenases, andslowing metabolism-related activities. Differential AS reinforced expression differences, with reducedAS frequencies typically observed for genes with upregulated expression. Significant differences werefound between the genotypes, including earlier and stronger upregulation of the alcohol fermentationpathway 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 offlooding and release from flooding stress contribute to increased flooding tolerance.Keywords: Helianthus annuus L.; crop; abiotic stress; flooding; alcohol dehydrogenase; transcrip-tomics; alternative splicing; co-expression network; mixed model; fermentation pathway1. IntroductionHistorically, people have inhabited regions with abundant water for multiple reasons,but most of all, for successful agriculture. However, the association with water means thatflooding typically is common as well, which can result in significant crop losses, especiallyfor crops that lack resistance to flooding. Model-based estimates of crop losses resultingfrom flooding stress suggest that such losses are likely to increase dramatically with climatechange 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 thelong-term goal of developing crops with greater resilience to flooding stress.From a plant physiological standpoint, flooding results in reduced gas exchange insubmerged roots and shoots, which can slow down rates of respiration and photosynthe-sis [3] or lead to damage at the cellular level (e.g., in lipids, proteins and DNA) due toincreased levels of reactive oxygen species (ROS) [4]. Plants that are frequently exposedAgronomy 2021, 11, 92. https://doi.org/10.3390/agronomy11010092 https://www.mdpi.com/journal/agronomyAgronomy 2021, 11, 92 2 of 21to flooding, however, have developed several strategies to respond and adapt to floodingstress, both physiologically and morphologically. These strategies include quiescence,anaerobic energy production, rapid shoot growth to escape flooding, and the developmentof aerenchyma to facilitate internal gas exchange.In the quiescence or “sit and wait” strategy, plants conserve energy such as ATP andcarbohydrates 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 withgrowth in 86 Arabidopsis accessions with natural variation for flooding tolerance, consistentwith the quiescence strategy [6].When oxygen, the essential substrate for respiration, becomes limited, plants are nolonger able to produce enough ATPs to meet cellular demands [7]. Under oxygen limitingconditions, regeneration of nicotinamide adenine dinucleotide (NAD) occurs via twofermentation pathways: (i) lactic acid fermentation and ii) alcohol fermentation [7–9]. Lacticacid fermentation results in the potentially damaging accumulation of lactate, whereasalcohol (ethanol) fermentation is considered to be less harmful. Therefore, the extent ofethanol production is thought to be positively correlated with a plant’s tolerance to floodingstress [10].Wide-ranging development of air-filled aerenchyma tissue can also contribute to flood-ing resistance by facilitating gas exchange between the shoot and the root [11]. However,this mechanism is ineffective for complete submergence, which typically triggers rapidshoot 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 rootscontaining aerenchyma represents another escape strategy by minimizing the distancebetween tissues that are in contact with air (above the waterline) and root tips [14–16].Adventitious root growth is promoted by increased ethylene gas [17,18] and is mediated bya combination of hormonal signals and reactive oxygen species [19]. The timing of adventi-tious root emergence is determined by the degree of flooding stress and the developmentstage of the plant, and it varies across plant species [16].Transcriptomic analysis can provide additional clues regarding the mechanisms under-lying resistance to flooding stress [20,21]. An example comes from soybean, which displaysextensive variation for waterlogging tolerance [22–26]. Transcriptional responses wereexplored by RNA-Seq, and flooding stress was associated with downregulation of photo-synthesis and chlorophyll synthesis related genes, consistent with reduced photosyntheticefficiency (i.e., quiescence) observed at the whole plant level [27]. Likewise, a combinationof transcriptome and proteome analyses showed that alcohol fermentation related geneswere significantly upregulated in soybean under flooding stress [28], indicating a role forfermentation as well. Arora et al. [29] conducted RNA-Seq in a flooding tolerant maizeinbred line under waterlogging and identified differentially expressed genes involved inseveral key pathways, including energy-production, programmed cell death, aerenchymaformation, and ethylene responsiveness, suggesting that multiple mechanisms underlieresponses to flooding stress. There is also evidence that some submergence-activated genesare clustered in the genome and that genes within these clusters are more highly expressedacross angiosperms than submergence-activated genes outside of such clusters [30].Alternative splicing (AS) offers a post-transcriptional mechanism for fine-tuningexpression 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, inanimals, exon skipping or shuffling is most frequent [32,33]. Most of these AS events(especially intron retention events) include premature stop codons, and transcripts areoften subjected to degradation via nonsense-mediated decay, leading to downregulation ofdosage [34–36]. Using RT–PCR, Syed et al. [37] found that in soybean, the SUB1 gene (amaster 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 appearsto contribute to homeostasis. In maize, inbred lines tolerant of waterlogging exhibiteddifferential splicing at prolyl 4-hydroxylases (which are thought to be oxygen sensorsAgronomy 2021, 11, 92 3 of 21under hypoxia stress), compared to intolerant lines when exposed to flooding stress [38].In Arabidopsis thaliana, transcriptomic analyses revealed water submergence-specific ASevents in pathways regulating energy reserves such as gluconeogenesis [39]. However, weare unaware of other studies that have investigated genome-wide AS under flooding stress.While transcriptomic responses to flooding stress have been investigated in a numberof crops, much remains to be understood about which responses are common acrossmany taxa and which are taxon-specific. In addition, previous studies typically havenot examined interactions of genotypes with environment and time to determine theidentity and timing of expression changes exclusive to tolerant lines. Here we reporton transcriptomic responses to flooding in sunflower (Helianthus annuus L.), which isone of the world’s most important oilseed crops, with global production valued at morethan 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, andflooding stress [40–44]. In North America, sunflowers sometimes experience short boutsof flooding in the spring, which results in the submergence of seeds, roots or even entireseedlings [45]. While such short-term flooding may not be lethal, yield and oil quality maybe reduced [40].Several potential adaptation strategies to flooding stress have been reported previ-ously for both wild and cultivated sunflowers, including aerenchyma development andanaerobic energy production. Wample and Reid [46] reported that flooding induces theproduction of adventitious roots de novo, as well as swelling (hypertrophy) of the lowerstem (hypocotyl) in cultivated sunflower. Both putative adaptations involve aerenchymaformation and are thought to represent a flooding escape strategy by promoting tissuegas exchange [46,47]. Likewise, evidence of increased anaerobic energy production hasbeen reported for wild sunflowers that occupy habitats that frequently experience springflooding [48]. A similar mechanism appears to be present in cultivated sunflower: floodingpromotes the synthesis of ethanol, which is successfully metabolized thanks to an increasein alcohol dehydrogenase (ADH) activity [49].Previously, we characterized phenotypic traits of 288 inbred sunflower lines from acultivated sunflower association mapping (SAM) population under both flooding stressand benign conditions [44]. We measured multiple traits for possible flooding responses,such as biomass and chlorosis reduction, hypocotyl hypertrophy and adventitious rootnumbers, and identified candidate genes underlying these responses. We found thatcultivated sunflower varied widely in flooding tolerance and that there was no inherenttradeoff between flooding tolerance and performance in the control treatment. On theother hand, we failed to find an effect of hypocotyl hypertrophy or adventitious rootdevelopment on flooding tolerance, perhaps because they represent a general responseto flooding across all cultivars. However, because of the ability of some of the genotypesto maintain growth without obvious negative impacts on photosynthesis (e.g., chlorosis),we suspected that increased anaerobic energy production might account for the success oftolerant phenotypes.In the current study, we have extended our analyses of flooding tolerance in sunflowerby studying expression and AS responses to flooding stress in one tolerant and susceptibleline. We identified the genes that are differentially expressed or differentially spliced inresponse 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 andpost-flooding. We also conducted a gene co-expression network analysis to identify ex-pression networks and to explore their connectedness and structure in resistant versussusceptible lines when responding to flooding stress. Our goal was not only to obtaina genome-wide assessment of molecular responses to flooding in sunflower but also toprovide additional clues regarding how flooding tolerance is achieved. We hypothesizethat pathways associated with aerenchyma formation and alcohol fermentation pathwaysplay key roles in general and specific responses to flooding, respectively.Agronomy 2021, 11, 92 4 of 212. Materials and Methods2.1. Plant Material and Flooding Stress TreatmentOut of 288 lines of the SAM population, we chose a resistant line (HA351) and asusceptible line (RHA428) based on vegetative growth performance and photosyntheticefficiency under the flooding stress (Figure 1). The plants were germinated and grownfollowing the protocol of Gao et al. [44]. Achenes were cleaned using 75% ethanol and0.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 thenplaced 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 theywere 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 weregrown 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 plantswere three weeks old, half were placed in large plastic bins filled with tap water to imposeflooding stress. The water level was kept at 1–2 cm above the soil surface throughout thestress treatment. After nine days, the pots were removed from the bins and returned topre-flooding conditions for six days to allow recovery. Leaf and root tissue were collectedfrom 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 dayof flooding), (ii) the last day (8th day of flooding) of the flooding treatment and (iii) 6 daysafter the flooding treatment ended (6th day of recovery). These time points were chosenbased on the results of Gao et al. [44].Agronomy 2021, 11, x FOR PEER REVIEW 5 of 22 (A) (B) Figure 1. Phenotypic comparisons between resistant (HA351) and susceptible (RHA428) lines. (A) Shoot biomass (g, dry weight) and (B) relative chlorophyll concentration of leaves (using a chlorophyll meter) from Gao et al. [44] that were performed across two cycles (CYI and CYII, independent experiments). 2.2. 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). Trimmo-matic (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 anal-ysis, HTSFilter (Version 1.12.0) [54] was used to remove genes with low constant expres-sion levels. Analyses of root and leaf samples were conducted separately due to the strong difference between the two (Figure S1). 2.3. 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 > log2FC or 1 < log2FC as calculated by edgeR (Version 3.18.1) [55], a Biocon-ductor 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 follow-ing 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 > log2FC or 1 < log2FC) 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 log2FC 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 pro-vided in Table S3. Figure 1. Phenotypic comparisons betwe n resistant (HA351) and susceptible (RHA428) lines. (A) Sho t biomas (g, dryweight) and (B) relative chlorophyll concentration of leaves (using a chlorophyll meter) from Gao et al. [44] that wereperformed across two cycles (CYI and CYII, independent experiments).A total of 96 sampl were collected: 2 genotypes × 2 tissues × 2 treatments × 3 timepoints × 4 biological replicates. For each biological replicate, leaflets and roots were pooledfrom three individual plants. Total RNA was extracted using the Ambion RNAqueouskit (AM1912, Invitrogen, Austin, TX, USA), and DNA contaminants were removed bythe 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 samplepreparation kit and sequenced on an Illumina HiSeq 4000 (PE75) by Genome Quebec inMontréal (Québec, Canada). Sequence data from this project have been deposited in theNational Center for Biotechnology Information Sequence Read Archive under accessionnumber PRJNA492303.A total of 48 seedlings (2 genotypes × 2 treatments × 3 time points × 4 biologicalreplicates) were measured following the protocols of Gao et al. [44]: shoot biomass (SB, gof dry weight), leaf biomass (LB, g of dry weight), leaf area (LA, cm2), leaf number (LN),Agronomy 2021, 11, 92 5 of 21plant height (PH, mm), stem diameter (SD, mm), chlorophyll concentration (CC, relativechlorophyll concentration unit via a chlorophyll meter), leaf mass per area (LMA, g/m2),hypocotyl diameter (HD, mm) and adventitious roots development (number of adven-titious roots, ARN) (Table S1). Raw trait values were then employed in a co-expressionanalysis (see below) to identify correlations with co-expression modules.2.2. RNA-Seq Data ProcessingThe quality of the 96 transcriptomes (75 bp paired-end raw reads, single-stranded) wasassessed 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]. Thefiltered reads were mapped to the sunflower reference genome [51] using STAR (Version2.5.2b) with non-default parameters of alignIntronMax 10,000 [52]. We mapped readsagainst the reference genome instead of a transcriptome assembly to avoid multi-mappingof reads onto different isoforms of the same gene, which is critical for accurate analysesof alternative splicing. Mapped reads were quantified by featureCounts in the Subreadpackage (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 expressionlevels. Analyses of root and leaf samples were conducted separately due to the strongdifference between the two (Figure S1).2.3. Analyses of Differential Gene Expression and Alternative SplicingFor each time point, we identified differentially expressed genes (DEG) betweenthe two conditions (control and flooding) using a false discovery rate (FDR) < 0.01 withfold-change of −1 > log2FC or 1 < log2FC as calculated by edgeR (Version 3.18.1) [55], aBioconductor package. Statistics of DEGs obtained from this approach are provided inTable S2. We converted raw read counts to reads per kilobase of transcript per millionmapped 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 fol-lowing the protocol of Tack et al. [57]. It identified the following major alternativesplicing events: intron retention (IR), an alternative donor (ALTD), alternative accep-tor (ALTA), alternative position (ALTP), exon skipping (SKIP), cryptic intron (CRIN) andcryptic exon (CREX). Reads mapping to the alternative and constitutive gene models werecounted 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 > log2FCor 1 < log2FC) to identify differentially spliced (DAS) genes since the proportion of alter-native splicing events is between 0 and 1 (calculated using a custom R script depositedon https://github.com/dejonggr/differential_as). We compared the log2FC between thetwo genotypes (pooling three time points together) to determine whether they differ intheir response to flooding stress. Statistics for DAS genes obtained from this approach areprovided in Table S3.2.4. Mixed Effect Model Analyses of Variation in Gene Expression and Alternative SplicingWe used generalized linear mixed effect models to identify genes that respondeddifferently in expression or AS to genotype, environment, time-point or the interactionbetween genotype and environment (G × E). We incorporated genotype, environment,time point and their interaction as fixed effects and individuals as random effects. Weused a gamma family distribution with an inverse link function and a Gaussian familydistribution with a logit transformation for expression (TPM value) and alternative splicingfrequency data (PSI value), respectively. We used an FDR to correct for multiple compar-isons 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 levelsbetween genotypes, environment or time points and corrected for multiple comparisonsusing the false discovery rate. All genes were tested for G × E in the mixed-effects model.Agronomy 2021, 11, 92 6 of 21In addition, we split each gene’s data into three time points. Each time point was tested forG × E in a fixed-effects model. All analyses were done using the programming languageR. Mixed effect models were performed using lme4 R package (Version 1.1-25) [60] andlmerTest (Version 3.1-3) [61] and analysis of variance was done using the “car” softwarepackage (Version 3.0-10) [62].2.5. Gene Co-Expression Network AnalysisWe 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 sunflowercultivars. The co-expression modules were assigned based on the TPM data of all expressedgenes 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. Phenotypictrait data (Table S1) were used to test for correlations with the transcriptome data. Theminimum 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 withthe highest positive absolute correlation coefficient value with phenotypic data, and thelowest P value, were selected for further analyses. Genes with high connectivity wereconsidered to be hub genes. Co-expression gene networks for the resistant and susceptiblelines were generated independently using the same parameters to permit comparison. Aconsensus network analysis was conducted using the same parameters as described aboveto identify consensus modules and to better assess overall similarities and dissimilaritiesbetween the resistant and susceptible line. Networks were illustrated using Cytoscape(Version 3.7.2) [64].2.6. Enrichment ComparisonsFor all gene sets, gene ontology (GO) enrichment analyses were performed usingagriGO Version 2.0 [65] with the Yekutieli multi-test adjustment method (false discoveryrate < 0.1) [66], a minimum of five mapping entries and Plant GO slim annotations. Theresults are shown via REVIGO [67] TreeMaps that summarize enriched GO terms and theirrelationships.3. Results3.1. RNA-Seq DataTo assess gene expression responses of cultivated sunflower to flooding stress, weconducted 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 thisexperiment (an average of 63 million reads per sample). After filtering out low-qualityreads and genes with low constant expression levels, 3,016,817,944 clean Illumina RNA-Seqreads for 48,611 expressed genes were used to identify DEGs.3.2. Differentially Expressed Genes (DEG) that Respond to Flooding Stress and Exhibit Genotype× Environment (G × E) InteractionsWe found 2338 genes (1510 upregulated, 828 downregulated) from leaf tissue and9089 genes (4417 upregulated, 4672 downregulated) from root tissue that were differentiallyexpressed under flooding stress (E) (FDR < 0.01) (Figure 2, Table S2). This root-dominantpattern of differential gene expression is shown in Table 2 and is displayed in Venn diagramsand scatter plots (Figure S2). In addition, more differentially expressed genes were foundduring the flooding periods (time-point 1 and 2) compared to the recovery period (time-point 3), especially in root tissue (Figure 2 and Figure S2, Table 2). Interestingly, thenumber 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 genesAgronomy 2021, 11, 92 7 of 21responding to the flooding stress and those at the recovery stage, may contribute to floodingtolerance (Table 2 and Table S2, Figure S2).Table 1. Different variables and comparisons employed in transcriptomic analyses.Variables ComparisonsTissue types Leaf vs. rootGenotypes (G) HA351 (resistant) vs.RHA428 (susceptible)Treatments (E) Control vs. floodingTime-points (1,2,3) 1st day of flooding vs. 8th day flooding vs. 6th day of recoveryG × E (HA351 vs. RHA428) × (control vs. flooding)G × E at 1,2 or 3 (HA351 vs. RHA428) × (control vs. flooding) at three differenttime-pointsExample: G × E1_rootThat is, DEGs from a mixed ANOVA of genotype (HA351 vs.RHA428) × treatment (control vs. flooding) interactions attime-point 1 in root tissue.Agronomy 2021, 11, x FOR PEER REVIEW 7 of 22 Table 1. Different variables and comparisons employed in transcriptomic analyses. Variables Comparisons Tissue types Leaf vs. root Genotypes (G) HA351 (resistant) vs. RHA428 (susceptible) Treatments (E) Control vs. flooding Time-points (1,2,3) 1st day of flooding vs. 8th day flooding vs. 6th day of recovery G × E (HA351 vs. RHA428) × (control vs. flooding) G × E at 1,2 or 3 (HA351 vs. RHA428) × (control vs. flooding) at three different time-points Example: G × E1_root That is, DEGs from a mixed ANOVA of genotype (HA351 vs. RHA428) × treatment (control vs. flood-ing) interactions a time- oint 1 i root tissue. 3.2. Differentially Expressed Genes (DEG) that Respond to Flooding Stress and Exhibit Geno-type × 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 differen-tially expressed under flooding stress (E) (FDR < 0.01) (Figu e 2, Table S2). This root-dom-inant 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 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 (Tables 2 and S2, Figure S2). Figure 2. A summary of differentially expressed genes that are responsive to flooding stress and exhibit interactions with genotype or time-point. E: treatment, G × E: genotype x treatment, G × E1,2 or 3: genotype x treatment at three different time-points. Figure 2. A summary of differentially expressed genes that are responsive to flooding stress andexhibit interactions with genotype or time-point. E: treatment, G × E: genotype × treatment,G × E1,2 or 3: genotype × treatment at three differe t time-points.Table 2. Number of differentially expressed genes (DEG) and differentially spliced genes (DAS) underflooding stress for the resistant (HA351) and susceptible (RHA428) cultivars at different time-points.Tissue Type Leaf RootCultivar: Time-PointDEG DAS DEG DASUp Down Up Down Up Down Up DownHA351: 1st day of flooding 6 35 22 136 781 554 220 136HA351: 8th day flooding 355 151 48 420 1178 1416 153 377HA351: 6th day of recovery 5 23 32 19 123 183 193 153RHA428: 1st day of flooding 22 48 11 2 1106 644 791 113RHA428: 8th day flooding 1006 294 31 183 2243 2825 23 233RHA428: 6th day of recovery 6 6 3 39 451 796 143 234Agronomy 2021, 11, 92 8 of 21Principle components analysis (PCA) offered further insights into gene expressionpatterns (Figure 3 and Figure S3). In root tissue, PC1 explained 23% of gene expressionvariation, which was partitioned mainly by treatment (control vs. flooding conditions).PC2 largely explained differences in genotype (resistant vs. susceptible lines), accountingfor 13.4% of the variation, whereas PC3 appears to distinguish time-point 1 from therest of the time-points and explains 10.4% of the variation. For leaf tissue, PC1 and PC2explained 30.1% and 13.5% of gene expression variation, respectively, but the variationwas not partitioned according to the fixed factors in our model. However, PC3 partitionssamples 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 contactof roots, but not leaves, to flooding stress. Indeed, the greater response of roots than leavesto flooding stress is seen throughout the experiment.Agronomy 2021, 11, x FOR PEER REVIEW 8 of 22 Table 2. Number of differentially expressed genes (DEG) and differentially spliced genes (DAS) under flooding stress for the resistant (HA351) and susceptible (RHA428) cultivars at different time-points. Tissue Type Leaf Root Cultivar: Time-Point DEG DAS DEG DAS Up Down Up Down Up Down Up Down HA351: 1st day of flooding 6 35 22 136 781 554 220 136 HA351: 8th day flooding 355 151 48 420 1178 1416 153 377 HA351: 6th day of recovery 5 23 32 19 123 183 193 153 RHA428: 1st day of flooding 22 48 11 2 1106 644 791 113 RHA428: 8th day flooding 1006 294 31 183 2243 2825 23 233 RHA428: 6th day of recovery 6 6 3 39 451 796 143 234 Principle components analysis (PCA) offered further insights into gene expression patterns (Figures 3 and S3). In root tissue, PC1 xplained 23% of gene expression variation, which was partitioned mainly by treatment (control vs. flooding conditio s). 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 be-tween 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. Figure 3. Principle component analysis (PCA) of differentially expressed genes under flooding stress in root tissue of sunflower cultivars. (A) PC1 vs. PC2. PC1 explains the difference between flooding stress and control/recovery as circles (control) and triangles (flooding) are clearly distinct, and PC2 explains the difference between genotypes as smaller shapes (resistant line) and larger shapes (susceptible line) are separated. (B) PC1 vs. PC3. PC1 explains the difference between flooding stress and control/recovery as circles (control) and triangles (flooding) are distinct, PC3 differentiates time-point 1 from the other time-points as red-colored shapes are distinguishable from green and blue colored ones, (C) Scree plot showing percentage variation explained by each PC. C = control, F = flooding; 1,2 or 3 refers to the different time points; HA351 = resistant line, RHA428 = susceptible line. 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 Figure 3. Principle compo ent analy is (PCA) of di ferentially s under flo ding stress in root tissue ofsunflower cultivars. (A) PC1 vs. PC2. PC1 explains the diff t een flooding stress and control/recovery as cir les(control) and triangles (flooding) are clearly distinct, and P ifference betwe n genotypes as smaller shapes(resistant line) and larger shapes (susceptible line) ar t . ( ) 1 vs. PC3. PC1 explains the differ nce betweenflooding stress and control/recovery as circles (control) and triangles (flooding) are distinct, PC3 differentiates time-point1 from the other time-points as red-colored shapes are distinguishable fro green and blue colored ones, (C) Scree plotshowing percentage variation explained by each PC. C = control, F = flooding; 1,2 or 3 refers to the different time points;HA351 = resistant line, RHA428 = susceptible line.To investigate which e ifferential expres ion between the resistant andsusceptible lines due to the fl t ent and/or the timing of the treatment, we useda mixed model approach (Figure 4; Table S2). After the FDR correction (0.1 or less), wefound that G× E interactions were significant for 1256 genes in leaf tissue and 1093 genes inroot tissue. Of these genes, we then searched for the subset that differed significantly (two-fold, FDR < 0.1) between lines under flooding stress but not under control conditions. Inleaf tissue (G × E_leaf), we found 8 genes that were significantly upregulated and 31 genesthat were significantly downregulated in the resistant line compared to the susceptibleline. In the root tissue (G × E root), we found 20 genes that were significantly upregulatedand 41 genes that were significantly downregulated in the resistant line relative to thesusceptible line. These genes represent good candidates for explaining the differencesbetween the two lines in their resistance to the flooding stress (Table S2–G × E leaf andG × E root).Agronomy 2021, 11, 92 9 of 21Agronomy 2021, 11, x FOR PEER REVIEW 9 of 22 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 (two-fold, 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 suscep-tible line. In the root tissue (G × E root), we found 20 genes that were significantly upreg-ulated 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). Figure 4. Example of results from mixed model analyses. (A) Heat map clustered by Euclidean distance showing top 10 differentially expressed genes (G × E1_root) with significantly higher expression in the resistant line compared to the susceptible line in root tissue. (B) HanXRQChr05g0161051 (alcohol dehydrogenase) was significantly upregulated under flooding stress at time-point 1 in root tissue of resistant line (G × E1_root). For the time-point analysis, we searched for the genes that were significantly differ-ent between lines at one or more time-points in the flooding treatment but not in the con-trol 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). 3.3. 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 co-dons, 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 deter-mine the frequency of reads that mapped to the alternatively spliced gene model com-pared 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 up-regulated, 760 downregulated) from leaf tissue and 2734 events (1578 upregulated, 1156 HanXRQChr05g0161051: alcohol dehydrogenase(B)(A)TPMLow HighControl Flooding Control Floodingtime-point1 time-point2 time-point3 time-point1 time-point2 time-point3time-point1 time-point2 time-point3 time-point1 time-point2 time-point31 2 3 1 2 3 1 2 3 1 2 3Row Z-score-4 -2 0 2 4Time-pointTreatmentSusceptible line Resistant lineGenotypeTop 10 most differentially expressed genes (G×E1_root)HanXRQChr13g0386591HanXRQChr05g0161051HanXRQChr06g0179461HanXRQChr17g0562911HanXRQChr06g0171031HanXRQChr11g0330061HanXRQChr13g0415181HanXRQChr04g0113431HanXRQChr10g0304671HanXRQChr08g0220261Resistant line Resistant lineSusceptible line Susceptible lineFigure 4. Example of results from mixed model analyses. (A) Heat map clustered by Euclidean distance showing top 10differentially expressed genes (G × E1_root) with significantly higher expression in the resistant line compared to thesusceptible line in root tissue. (B) HanXRQChr05g0161051 (alcohol dehydrogenase) was significantly upregulated underflooding stress at time-point 1 in root tissue of resistant line (G × E1_root).For the time-point analysis, we searched for the genes that were significantly dif-ferent between lines at one or more time-points in the flooding treatment but not in thecontrol treatment (Table S2). In leaf tissue, we found 91 genes that were upregulated and178 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 susceptibleline. 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 theresistant versus susceptible line (G × E1-3_root). The identity and function of these genesmay offer clues regarding how shifts in the timing of expression of particular categories ofgenes contribute to flooding resistance (Table S2– G × E1-3_leaf and G × E1-3_root).3.3. Differentially Spliced Genes (DAS) That Respond to the Flooding StressWe identified AS events and estimated their frequencies using the transcriptome datadescribed above (Table S4). As xpected, intron retention events w re most frequ nt andexon 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) anddownre ulate the overall expression of the ge e. An AS ratio was calculated to determinethe frequ ncy of reads th t ma ped to the alternatively spliced gene model compared toth constitutive gene model. Thus, upregulation o downregulation of AS ref rs to anincrease or decrease, r sp ctively, of the AS ratio. We f und 892 AS vents (132 upregulated,760 downr gulated) from leaf tissue and 2734 events (1578 pregulated, 1156 downregu-lated) from root tissue th varied significantly in frequency u der flooding stress (Table S3).Similar to the DEGs described above, differentially spliced events were more than twice ascommon in root compared to leaf tissue. A total of 69 genes (11 upregulated, 58 downregu-lated) from leaf tissue and 97 genes (56 upregulated, 41 downregulated) from root tissuewere 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 otherhand, PC2 of leaf tissue (7.6% of variation) and PC3 of root tissue (8.4% of variation) wereassociated with the genotypic differences.3.4. Putative Functions of Differentially Expressed GenesOf the genes that were differentially expressed under flooding stress (E leaf), thoseupregulated in leaf tissue were enriched in 71 GO terms (Table S5). Many of the significantAgronomy 2021, 11, 92 10 of 21GO terms were related to protein translation and response to an extracellular stimulus(biological process; Figure S5A), ribosome (cellular component, Figure S5B), and structuralmolecule activity (molecular function; Figure S5C). Fifty-five GO terms were significantlyenriched for genes that were downregulated in leaf tissue (E_leaf); GO terms relatingto secondary metabolism, homeostatic process and photosynthesis (biological process;Figure S5D), cell wall and membrane (cellular component, Figure S5E), and catalytic activ-ity and binding (molecular function; Figure S5F) were especially frequent. Overall, theseresults indicate that in leaf tissue under flooding stress, genes involved in the productionof extracellular stimulus-responsive peptides (possibly due to ethylene and ROS) wereupregulated, whereas those related to secondary metabolism, energy production or photo-synthesis were downregulated. Slowing down these and other metabolic pathways maycontribute to flooding tolerance.For root tissue under flooding stress, upregulated genes (E_root) were enriched in75 GO terms (Table S5). The most significant GO terms were enriched in response to anextracellular stimulus (biological process; Figure S5G), cell wall and membrane (cellularcomponent, Figure S5H), catalytic activity and transferase activity (molecular function;Figure S5I). Genes that were significantly downregulated in root tissue (E_root) wereenriched 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, roottissue responded to flooding stress as expected, by upregulating expression levels offlooding responsive genes. In addition, many metabolism-related activities appear to haveslowed, presumably due to a lack of oxygen and energy.We also conducted GO enrichment analyses for the DEGs identified in the mixedmodels (Table S5). As a reminder, in these analyses, upregulated genes are those that areupregulated in the resistant line under flooding conditions when compared to the suscep-tible line in either treatment and to the resistant line in the control treatment. Likewise,downregulated genes are those that are downregulated in the resistant line in the floodingtreatment relative to the susceptible line in either treatment and to the resistant line undercontrol conditions.In the leaf tissue (G × E_leaf), only eight genes with significant G × E were upregu-lated, which is too few for GO enrichment analysis. Nonetheless, several of the genes werenotable and include HanXRQChr08g0214611 (zinc finger, RING/FYVE/PHD-type knownas XERICO) and HanXRQChr14g0450491 (high-affinity nitrate transporter 2.7). Thirty-onegenes with significant G × E interactions were downregulated; this set of genes was en-riched in response to stress (biological process), extracellular region (cellular component),and catalytic activity and binding (molecular function). Results of G × E analyses weresimilar for each time point (G × E1-3_leaf): signaling-related genes were upregulated, andresponse to stress genes were downregulated. This implies that stress response genes inthe resistant line are less responsive to flooding than in the susceptible line, which mayunderlie the greater growth performance of the former.In the root tissue (G × E_root), upregulated genes exhibiting significant G × E interac-tions were related to GO: 0003824 “catalytic activity” (FDR = 0.031; molecular function) andGO: 0005576 “extracellular region” (FDR = 0.049; cellular component), whereas, downregu-lated 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 down-regulated 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 andGO: 0003824 “catalytic activity” (FDR = 8.00E-08) in the downregulated fraction, andsignificant genes at time-point 3 (G × E3_root) were enriched in GO:0050896 “responseto stimulus” (FDR = 3.60E-05) in the upregulated subset and GO:0003824 “catalytic ac-tivity” (FDR = 0.001) in the downregulated subset. These enrichment patterns indicatethat the resistant and susceptible line differ significantly both in the kinds of genes thatAgronomy 2021, 11, 92 11 of 21respond to flooding stress and the timing of the response. Genes involved in catalyticactivity appear to play an especially critical role in flooding response. Notably, four ADHgenes (HanXRQChr05g0161031, HanXRQChr05g0161051, HanXRQChr06g0179461 andHanXRQChr09g0243891) were upregulated at time-point 1 (G × E1_root), as well as twoadditional ADH genes (HanXRQChr05g0161051 and HanXRQChr17g0559391) at time-point 2 (G × E2_root). This implies that rapid upregulation of the fermentation pathwaycontributes to the greater flooding tolerance of the resistant line.3.5. Putative Function of Differentially Spliced GenesAmong 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 carbohydratemetabolism (biological process; Figure S6A), intracellular organelle part (cellular compo-nent, Figure S6B), and pyrophosphatase activity (molecular function; Figure S6C). Geneswith decreased AS in leaf tissue were significantly enriched in embryo development (bio-logical process; Figure S6D), intracellular organelle part (cellular component, Figure S6E),and in transferase activity (molecular function; Figure S6F). Increased AS frequency mayreduce the expression level of metabolic genes via NMD, further downregulating metabolicpathways 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 (cellularcomponent, Figure S6H), and nuclease activity (molecular function; Figure S6I). Genes withdecreases 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 embryodevelopment are accord with those in leaf tissue, suggesting the regulation of AS in bothleaf and root tissue work in a similar way and likely acts to fine-tune gene expressionresponses to environmental change.The mixed model analysis (Table S3) indicated that very few AS events exhibit sig-nificant G × E interactions, so we were unable to be able to conduct GO enrichmentanalyses. Therefore, individual DAS events were explored. For example, AS frequenciesof HanXRQChr09g0273081 (elongation factor 1-delta) and HanXRQChr11g0325531 (trans-lation 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 roottissue, HanXRQChr05g0161051 (alcohol dehydrogenase) had reduced AS in the resistantline at time-point 1 (G × E1_root), and this may contribute to more efficient expressionvia constitutive splicing in the fermentation pathway. HanXRQChr12g0371921 (myc-type, basic helix-loop-helix) showed an increase in AS in the resistant line at time-point1 (G × E1_root), and its molecular function (GO: 0003700) is related to DNA-bindingtranscription factor activity.Overall, the changes in AS reinforce expression changes under flooding stress, withstrong negative correlations observed (r ranges from−0.57 to−0.61 and P from 5.01 × 10−11to 1.41 × 10−49) regardless of tissue type or treatment (Figure S11).3.6. Gene Co-Expression Network Modules Underlying Responses to Flooding StressCo-expression networks were constructed using all 96 RNA-Seq transcriptomes, butroot and leaf transcripts were analyzed separately due to the large differences betweenthe two tissue types. We generated network modules for resistant and susceptible linesindependently using standard parameters (see Materials and Methods) and comparedthem via consensus analysis following the protocol of Langfelder and Horvath [63]. Wewere particularly interested in correlations that differed significantly between resistant andsusceptible data sets in strength or sign, thereby permitting identification of unique geneco-expression networks in the flooding resistant sunflower cultivar.The leaf tissue network for the resistant line was comprised of 34 clusters of co-expressed genes. These module eigengenes (ME) were identified by the first componentAgronomy 2021, 11, 92 12 of 21in a PCA, which was generated by the WGCNA package. Module-trait associations arepresented by a heat map (Figure S7A), which illustrates the correlations between theMEs and phenotypic traits. Chlorophyll concentration is closely connected to floodingresistance in sunflower due to its impact on photosynthesis [44]. We observed that genesin module ME-black (Table S7) in the heat map are most positively correlated (correlationcoefficient = 0.57, P = 0.004) with chlorophyll concentration. ME-black consists of genes thatwere 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 theconsensus data analysis, which shows substantial overlap with the ME-black module fromthe 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 anon-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 contributeto 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 differ-entially expressed genes identified by the mixed model, and if so, how they are connected orembedded 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 differ-entially expressed genes in the G × E_leaf analysis (Table S7; Figure S7E). Many of the hubgenes 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 down-regulated) 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 3-related 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 genecontinued 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 downregu-lated in the resistant versus susceptible line (Table S2), consistent with our hypothesis thatthe resistant line recovers more quickly than the susceptible line by turning down stressresponse 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 nothave root phenotypic trait data, any correlations between root tissue expression and phe-notypic changes are expected to be indirect (e.g., development or growth rate). We chose tofocus on the ME-red module (Table S7) since it has the strongest positive correlation amongmodule-trait relationships, and it is positively correlated with leaf biomass (correlationcoefficient = 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-redconsists of the genes that were enriched in GO terms (Table S8) related to the response to en-dogenous 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 (correla-tion coefficient = 0.55, P = 0.005) with the leaf area in the resistant line (Figure S10B), butnot in the susceptible line (correlation coefficient = −0.04, P = 0.9; Figure S1C). The lackof consensus for this trait between resistant and susceptible lines is indicated as NA inFigure S10D.The ME-red module in the root tissue of the resistant line was further analyzed byidentification of hub genes. As before, we used Cytoscape for network visualizationwith the top 10 hub genes (with the highest degree of connection), as well as differen-Agronomy 2021, 11, 92 13 of 21tially expressed genes from the G × E analysis (Table S7; Figure 5). The hub genes withthe highest degree of connectivity include HanXRQChr09g0248141 (hypoxia-responsivefamily protein), HanXRQChr12g0380481 (transcription factor TGA like domain), andHanXRQChr15g0464051 (Thiamine pyrophosphate dependent pyruvate decarboxylasefamily protein). Twenty-one differentially expressed genes exhibiting significant G × Einteractions 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 describedpreviously, suggesting that more rapid and stronger upregulation of the fermentationpathway in the resistant line may contribute to its greater tolerance to flooding.Agronomy 2021, 11, x FOR PEER REVIEW 13 of 22 among module-trait relationships, and it is positively correlated with leaf biomass (corre-lation 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 correla-tion (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 ex-pressed genes from the G × E analysis (Table S7; Figure 5). The hub genes with the highest degree of connectivity include HanXRQChr09g0248141 (hypoxia-responsive family pro-tein), 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 A H homologs (HanXRQChr05g0161031, HanXRQChr05g0161051, HanXRQChr06g0179461) were found. One of these (HanXRQChr05g0161051) was also upreg lated at time-point 2 (G × E2_root) (Table S2). This is consistent with the results from the ixed odel analyses de-scribed previously, s ggesting that more rapid and stronger upregulation of the fermen-tation pathway in the resis ant line may contribute to its greater tol rance to flooding. Figure 5. Results of gene co-expression network analysis in root tissue. (A) Venn diagram of ME-red and differentially expressed genes with significant two-way genotype by treatment (i.e., envi-ronment) effects. GE indicates genotype x 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 ex-pressed genes with significant genotype by treatment interactions when combined across the three time points (except for the ADH homologs). Figure 5. Results of gene co-expression network analysis in root tissue. (A) Venn diagram of ME-red and differ-entially expressed genes with significant two-way genotype by treatment (i.e., environment) effects. GE indicatesgenotype × environment interactions, and the number indicates the time-point. (B) Visualization of network connec-tions among the most connected genes in ME-red, as generated by Cytoscape software. Orange circles are hub geneswith the highest connectivity. Green circles are alcohol dehydrogenase (ADH) homologs. White circles are differentiallyexpressed genes with significant genotype by treatment interactions when combined across the three time points (except forthe ADH homologs).3.7. Comparison with Evolutionarily Conserved Submergence Activated GenesA comprehensive study of submergence upregulated gene families (SURFs) wasreported recently by Reynoso et al. [30]. They analyzed root tissue of four flowering plantspecies: rice, legume (Medicago truncatula), domesticated tomato (Solanum lycopersicum)and its dryland-adapted wild relative (Solanum pennellii) to understand variation andconservation of SURFs among species.In our root expression data, we identified 35 homologs of conserved SURFs that re-sponded significantly to flooding, of which 28 were upregulated and seven downregulated(Table S2). Two of these SURFs (AT1G77120—alcohol dehydrogenase and AT4G33070—pyruvate decarboxylase 1) were found to upregulated in the resistant line relative to thesusceptible line at time-points 1 and 2, and at time-point 1, respectively (Table S2). Pyruvatedecarboxylase initiates the ethanol fermentation process under anaerobic conditions byconverting pyruvate into acetaldehyde and carbon dioxide, providing additional evidencefor the importance of accelerated upregulation of the alcohol fermentation pathway.Agronomy 2021, 11, 92 14 of 21With respect to alternative splicing, a significant increase in AS frequency was foundfor six homologs of SURF genes, whereas for three genes, a significant decrease in ASfrequency was observed (Table S3). Notably, one of the genes with reduced AS wasHanXRQChr05g0161021 (AT1G77120—alcohol dehydrogenase), further increasing thenumber of functional transcripts of this gene under flooding stress.3.8. Case Study of the Alcohol Fermentation PathwayGiven previous evidence pointing towards the importance of the alcohol fermentationpathway in flooding tolerance in sunflower [49], as well as differences between the resistantand susceptible line observed in the present study, we decided to take a more comprehen-sive 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 (mostlyupregulated) under flooding stress; however, none had significant G × E interactions(Table S9), suggesting that upregulation of ADH homologs is a general response to floodingstress in leaf tissue. In contrast, many ADH homologs were differentially expressed underflooding stress (both up and down) in the root tissue (Table S9), including five of whichwere upregulated in the resistant line, but not the susceptible line. Four of these were up-regulated at time-point 1 (G × E1_root) (HanXRQChr05g0161031, HanXRQChr05g0161051,HanXRQChr06g0179461, and HanXRQChr09g0243891) (Table S9), and three were includedin the red module of the resistant line that was described in the co-expression networkanalysis.Noticeably, four ADH homologs that were responsive to flooding appear to betandemly duplicated on chromosome 5 (Figure 6A). All four had higher expression in theresistant versus susceptible line under the flooding stress, with the strongest differenceobserved 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 linerelative to the sensitive line at time-point 1 (G × E1_root; FDR = 9.95 × 10−6; Figure 6B).Agronomy 2021, 11, x FOR PEER REVIEW 15 of 22 Figure 6. Alcohol dehydrogenase (ADH) homologs that are differentially expressed or spliced. (A) ADH homologs on Chromosome 5 in root tissue that were upregulated by flooding stress. (B) The frequencies of an alternative splicing event (alternative acceptor) between exon 1 and exon 2 of an ADH homolog (HanXRQChr05g0161051) displayed significant genotype × environment interactions. The susceptible line exhibited greater AS in the flooding treatment, resulting in fewer functional transcripts compared to the resistant line. 4. Discussion 4.1. 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 sev-eral other taxa, including Taxodium, a flooding-tolerant conifer. Taxodium plantlets exhib-ited 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 shoot-dominant ex-pression 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. 4.2. General Responses to Flooding We observed two main responses to flooding stress (E_leaf and root) that were com-mon to both the resistant and susceptible lines: upregulation of ethylene-responsive tran-scription factors and downregulation of metabolic, energy production, and photosynthe-sis-related pathways. The accumulation of ethylene gas is considered to be the most con-sistent, timely and reliable indicator of early flooding stress [70], and ethylene-responsive transcription factors (ERFs) play an early and critical role in mediating transcriptomic re-sponses to flooding stress [71]. Ethylene regulates a hormonal cascade (abscisic acid, gib-berellic 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 signifi-cantly 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 i re 6. lc l e r e ase ( ) l s t at are differentiall expresse or spliced. ( ) ho l s oChromosome 5 in root tissue that were upregulated by flooding stress. (B) The frequencies of an alternative splicing event(alternative acceptor) between exon 1 and exon 2 of an ADH homolog (HanXRQChr05g0161051) displayed significantgenotype × environment interactions. The susceptible line exhibited greater AS in the flooding treatment, resulting in fewerfunctional transcripts compared to the resistant line.. Discussion4.1. Root Dominant Responses to Flooding Stress in SunflowerChanges in expression and AS were greater in r ot tha shoot tissue in all analysesand at all time-points. This root-dominant response to flooding has been observed inseveral other taxa, includi g Taxodium, a fl oding-tolerant conif r. Taxodium plantletsexhibited root dominant differential expression of the antioxidative defense, glycolysis, a dAgronomy 2021, 11, 92 15 of 21fermentation pathways, which play a key role in adaptive flooding responses by reducingdamage from ROS species in the former pathway and maintaining ATP production inthe latter two pathways [68]. Likewise, gray poplar (Populus × canescens) displayedgreater transcriptional change in roots than leaves under flooding stress and upregulatedgene pathways were similar to those reported for Taxodium [69]. However, there areexceptions to this pattern. For example, van Veen et al. [39] reported greater expressiondivergence in shoot than root tissue in Arabidopsis under flooding stress. This may bea consequence of different resistance strategies, with Arabidopsis exhibiting rapid shootgrowth to escape submergence, whereas sunflower (for example) appears to upregulateanaerobic energy production in roots and produces adventitious roots. In contrast to shoot-dominant expression divergence, Arabidopsis shows greater differential AS in roots whenexposed to flooding stress [39], an observation that aligns with our study. Establishing rootdominant AS as a general rule, however, awaits other studies of AS responses to floodingacross multiple tissue types.4.2. General Responses to FloodingWe observed two main responses to flooding stress (E_leaf and root) that were com-mon to both the resistant and susceptible lines: upregulation of ethylene-responsive tran-scription factors and downregulation of metabolic, energy production, and photosynthesis-related 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 transcrip-tion factors (ERFs) play an early and critical role in mediating transcriptomic responsesto flooding stress [71]. Ethylene regulates a hormonal cascade (abscisic acid, gibberellicacid and auxin) that along with ROS promote adventitious root growth [11], and severalprevious studies have shown that upregulation of ERFs is positively correlated with ad-ventitious root formation (e.g., in cucumber [72] and rice [73]). ERFs were significantlyupregulated 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 stressin our phenotypic data (Table S1), as well as in our previous large-scale GWA study ofsunflower cultivars [44]. It is important to keep in mind, though, that upregulation of ERFsand increased adventitious rooting are seen in both the resistant and susceptible lines, sothis 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 tolower availability of oxygen, which can lead to chlorosis, leaf abscission and a declinein metabolic rates [12,20,74,75]. Better control of metabolic pathways and energy usecan reduce negative impacts from flooding stress [15,76]. In sunflower, we observedphenotypic evidence of reduced photosynthesis activity (Figure 1, Table S1), and similarresults were reported previously by Gao et al. [44]. Reduced photosynthetic activity can alsobe inferred from the downregulation of metabolic, energy production, and photosynthesis-related pathways in both the resistant and susceptible cultivars in the flooding treatment(Figure S5D). The downregulation of photosynthesis-related gene pathways followingflooding has previously been reported in several other plant species, including soybean [27]and canola [77].4.3. Changes in AS Levels Reinforce Expression DifferencesAlternative splicing is a post-transcriptional mechanism that can increase the diversityof the proteome and regulate gene expression [78–80]. The latter typically occurs via theexpression of alternative transcripts that are targeted to NMD [81,82]. Increases in thefrequency 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 nonsenseisoforms (and a lower proportion of functional transcripts) in Arabidopsis thaliana followingtemperature stress. On the other hand, decreasing frequencies of AS can restore thefunction of genes that are inhibited by AS. For example, the DEHYDRATION-RESPONSIVEAgronomy 2021, 11, 92 16 of 21ELEMENT BINDING 2 (DREB2) transcription factor produces a non-functional transcriptunder 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 werestrongly negatively correlated with levels of gene expression (Figure S11), thereby rein-forcing expression divergence. For example, metabolism-related gene pathways exhibiteddecreased 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 hadsignificant differences in both expression and AS was relatively small, in most cases, dif-ferences were in the opposite direction, leading to reinforcing effects on the number offunctional transcripts. This was most pronounced in leaf tissue (E_leaf), where genes thatshowed increased expression and downregulated AS were most frequent (Figure S12).Such a complementary pattern between expression and AS levels has been previouslyreported in a study of homologous expression and AS patterns in canola under severalabiotic stresses [86]. A focus of future studies should be to validate these inferences usingproteomic approaches.4.4. Expression Changes Specific to the Flooding Resistant CultivarThe flooding resistant cultivar differed from the susceptible cultivar in three mainways: earlier and stronger upregulation of the alcohol fermentation pathway mutedexpression of stress-responsive genes in leaf tissue and more rapid recovery of pre-stressexpression levels. Below we discuss how these differences may lead to greater resistanceor tolerance.The alcohol fermentation pathway represents a critical metabolic, and physiologicalresponse mechanism to hypoxia in higher plants by providing the energy (NAD+) neededfor glycolysis to continue in the absence of oxygen [7–9,87]. The pathway is a highlyconserved adaptation to low oxygen, and upregulation of the pathway can enhance aplant’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 bothdetoxifies toxic acetaldehyde and regenerates oxidized NAD+ [87]. We explored the ex-pression of ADH homologs (orthologous to AT1G77120 in Arabidopsis thaliana). Our mixedmodel analysis of genotype x treatment across the three different time points revealed thata number of ADH homologs showed earlier and/or stronger expression in the resistantthan susceptible genotype when under flooding stress (G × E1_root vs. G × E2_root orG × E3_root). These differences were further reinforced by significant downregulationof ADH AS levels (leading to a larger fraction of functional transcripts) in the resistantline relative to the susceptible line. Notably, five ADH homologs on chromosome 5 ap-pear to be the product of tandem duplication, and they were upregulated much morestrongly 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 tocontribute to abiotic stress tolerance in plants. A well-known example is the polygenicSUBMERGENCE 1 (SUB1) locus, which confers tolerance to flooding in rice [90]. In ad-dition, most ADH homologs were found in a flooding resistance-related co-expressionnetwork in root tissue, further supporting the importance of this gene family. The genesfound in this co-expression network largely overlap with hypoxia resistance networkspreviously reported for other plants [39,91], except for the inclusion of alcohol fermentationpathway in the sunflower network. The hub genes in this co-expression network, such asHanXRQChr08g0216771 (phosphofructokinase) and HanXRQChr12g0380481 (transcriptionfactor TGA), likely represent key regulators of the alcohol fermentation pathway, althoughthis 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 ofthe resistant genotype and that the resistant genotype recovers pre-stress expression levelsmore quickly than the sensitive genotype. The rapid recovery of metabolic homeostasis isAgronomy 2021, 11, 92 17 of 21often associated with stress tolerance [92,93], but whether this is a cause or consequence ofincreased resistance is not clear in most cases. In sunflower, it may be that the sensitivegenotype over-responds to flooding, leading to reduced growth rates. This has beenreported for rice, in which overexpression of a stress-responsive gene, OsHOX24, reducesroot and shoot growth [94]. Alternatively, the greater energy production in roots of theresistant genotype (and/or other changes not discussed here) may reduce the negativeimpacts of flooding, leading to a less severe response and more rapid recovery.5. ConclusionsWe conducted a large-scale transcriptomics study to investigate the eco-physiologicalresponses of sunflower, an important oilseed crop, to flooding stress. Using mixed modeland co-expression network approaches, we explored the differential expression and ASpatterns to identify genes and genetic pathways that underlie general responses to flood-ing stress, as well as those specific to a resistant cultivar. We found that upregulationof ethylene-responsive transcription factors and adventitious root growth occurred inboth the sensitive and resistant cultivated lines and appear to represent a general criticalresponse to flooding stress. Our hypothesis that differential expression and AS of thealcohol fermentation pathway contributes to the greater tolerance of the resistant linewas supported. Greater flooding tolerance of the resistant line is also associated withmuted expression of stress-responsive genes in leaf tissue and with a more rapid recoveryof “normal” expression after the stress period is over. A goal for the future will be toextend our results across more genotypes and to validate expression and AS patterns withproteomic data.Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4395/11/1/92/s1, Figure S1: Venn diagrams of differentially expressed or spliced genes underflooding stress in leaf versus root tissue; Figure S2: Scatter plot of the log2 ratio versus meancount per million (CPM); Figure S3: Principle component analysis (PCA) of differentially expressedgenes under flooding stress in leaf tissue of sunflower cultivars; Figure S4: Principle componentanalysis (PCA) of differentially spliced genes under flooding stress in leaf tissue (A-C) and roottissue (D-F) of sunflower cultivars; Figure S5: Gene ontology (GO) enrichment analysis resultsof expression data represented by the TreeMap view of REVIGO; Figure S6: Gene ontology (GO)enrichment analysis results of alternative splicing data represented by the TreeMap view of REVIGO;Figure S7: Gene co-expression network analysis of leaf tissue of resistant line; Figure S8: Consensusmodules in leaf tissue; Figure S9: Gene co-expression network analysis of root tissue of resistantline; Figure S10: Consensus modules in root tissue; Figure S11: Scatter plot of Pearson’s productmoment 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: Differentiallyexpressed gene responding to different variables; Table S3: Differentially spliced gene responding todifferent variables; Table S4: Frequency of alternative splicing events estimated from the sunflowertranscriptome data reported here; Table S5: Gene ontology enrichment analysis results of differentiallyexpressed 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 resultsof gene co-expression network modules; Table S9. Alcohol dehydrogenase homologs in sunflower.Author Contributions: Conceptualization, J.S.L., L.G. and L.H.R.; methodology, J.S.L., L.G., L.M.G.and L.H.R.; software, J.S.L. and L.M.G.; data curation, J.S.L.; writing—original draft preparation,J.S.L., and L.H.R.; writing—review and editing, J.S.L., L.G., L.M.G. and L.H.R.; visualization, J.S.L.;supervision, L.H.R.; project administration, J.S.L. and L.H.R.; funding acquisition, L.H.R. All authorshave read and agreed to the published version of the manuscript.Funding: Funding was provided by Genome BC and Genome Canada (LSARP2014-223SUN).Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.Agronomy 2021, 11, 92 18 of 21Data Availability Statement: Data presented in this study are openly available in the National Centerfor Biotechnology Information Sequence Read Archive under accession number PRJNA492303.Conflicts of Interest: The authors declare no conflict of interest.References1. Rosenzweig, C.; Tubiello, F.N.; Goldberg, R.; Mills, E.; Bloomfield, J. Increased crop damage in the US from excess precipitationunder climate change. Glob. Environ. Chang. 2002, 12, 197–202. [CrossRef]2. Kaniewski, D.; Marriner, N.; Morhange, C.; Faivre, S.; Otto, T.; Van Campo, E. Solar pacing of storm surges, coastal flooding andagricultural losses in the Central Mediterranean. Sci. Rep. 2016, 6, 25197. [CrossRef] [PubMed]3. Bailey-Serres, J.; Voesenek, L.A.C.J. Flooding stress: Acclimations and genetic diversity. Annu. Rev. Plant Biol. 2008, 59, 313–339.[CrossRef] [PubMed]4. Fukao, T.; Yeung, E.; Bailey-Serres, J. The submergence tolerance regulator SUB1A mediates crosstalk between submergence anddrought tolerance in rice. Plant Cell 2011, tpc-110. [CrossRef] [PubMed]5. Colmer, T.D.; Voesenek, L.A.C.J. Flooding tolerance: Suites of plant traits in variable environments. Funct. Plant Biol.2009, 36, 665–681. [CrossRef]6. Vashisht, D.; Hesselink, A.; Pierik, R.; Ammerlaan, J.M.H.; Bailey-Serres, J.; Visser, E.J.W.; Pedersen, O.; van Zanten, M.;Vreugdenhil, D.; Jamar, D.C.L.; et al. Natural variation of submergence tolerance among Arabidopsis thaliana accessions. N. Phytol.2011, 190, 299–310. [CrossRef]7. Van Dongen, J.T.; Licausi, F. Oxygen sensing and signaling. Annu. Rev. Plant Biol. 2015, 66, 345–367. [CrossRef]8. Tadege, M.; Dupuis, I.; Kuhlemeier, C. Ethanolic fermentation: New functions for an old pathway. Trends Plant Sci.1999, 4, 320–325. [CrossRef]9. Yang, W.; Catalanotti, C.; Posewitz, M.C.; Alric, J.; Grossman, A.R. Insights into algal fermentation. In Low-Oxygen Stress in Plants;van Dongen, J., Licausi, F., Eds.; Springer: Vienna, Austria, 2014; Volume 21, pp. 135–163. [CrossRef]10. Blokhina, O.; Virolainen, E.; Fagerstedt, K.V. Antioxidants, oxidative damage and oxygen deprivation stress: A review. Ann. Bot.2003, 91, 179–194. [CrossRef]11. Loreti, E.; van Veen, H.; Perata, P. Plant responses to flooding stress. Curr. Opin. Plant Biol. 2016, 33, 64–71. [CrossRef]12. Voesenek, L.A.C.J.; Bailey-Serres, J. Flood adaptive traits and processes: An overview. N. Phytol. 2015, 206, 57–73. [CrossRef][PubMed]13. Bailey-Serres, J.; Voesenek, L.A.C.J. Life in the balance: A signaling network controlling survival of flooding. Curr. Opin. PlantBiol. 2010, 13, 489–494. [CrossRef] [PubMed]14. Voesenek, L.A.C.J.; Bailey-Serres, J. Flooding tolerance: O2 sensing and survival strategies. Curr. Opin. Plant Biol. 2013, 16, 647–653.[CrossRef] [PubMed]15. Bailey-Serres, J.; Lee, S.C.; Brinton, E. Waterproofing crops: Effective flooding survival strategies. Plant Physiol. 2012, 160, 1698–1709.[CrossRef] [PubMed]16. Steffens, B.; Rasmussen, A. The physiology of adventitious roots. Plant Physiol. 2016, 170, 603–617. [CrossRef]17. Kende, H.; Van Der Knaap, E.; Cho, H.T. Deepwater rice: A model plant to study stem elongation. Plant Physiol.1998, 118, 1105–1110. [CrossRef]18. Lorbiecke, R.; Sauter, M. Adventitious root growth and cell-cycle induction in deepwater rice. Plant Physiol. 1999, 119, 21–30.[CrossRef]19. Steffens, B.; Kovalev, A.; Gorb, S.N.; Sauter, M. Emerging roots alter epidermal cell fate through mechanical and reactive oxygenspecies signaling. Plant Cell 2012, 24, 3296–3306. [CrossRef]20. Bailey-Serres, J.; Fukao, T.; Gibbs, D.J.; Holdsworth, M.J.; Lee, S.C.; Licausi, F.; Perata, P.; Voesenek, L.A.C.J.; van Dongen, J.T.Making sense of low oxygen sensing. Trends Plant Sci. 2012, 17, 129–138. [CrossRef]21. Mustroph, A. Improving flooding tolerance of crop plants. Agronomy 2018, 8, 160. [CrossRef]22. VanToai, T.T.; Beuerlein, A.F.; Schmitthenner, S.K.; St. Martin, S.K. Genetic variability for flooding tolerance in soybeans. Crop. Sci.1994, 34, 1112–1115. [CrossRef]23. VanToai, T.T.; St. Martin, S.K.; Chase, K.; Boru, G.; Schnipke, V.; Schmitthenner, A.F.; Lark, K.G. Identification of a QTL associatedwith tolerance of soybean to soil waterlogging. Crop. Sci. 2001, 41, 1247–1252. [CrossRef]24. Reyna, N.; Cornelious, B.; Shannon, J.G.; Sneller, C.H. Evaluation of a QTL for waterlogging tolerance in southern soybeangermplasm. Crop. Sci. 2003, 43, 2077–2082. [CrossRef]25. Valliyodan, B.; Ye, H.; Song, L.; Murphy, M.; Shannon, J.G.; Nguyen, H.T. Genetic diversity and genomic strategies for improvingdrought and waterlogging tolerance in soybeans. J. Exp. Bot. 2016, 68, 1835–1849. [CrossRef] [PubMed]26. Wu, C.; Zeng, A.; Chen, P.; Florez-Palacios, L.; Hummer, W.; Mokua, J.; Klepadlo, M.; Yan, L.; Ma, Q.; Cheng, Y. An effective fieldscreening method for flood tolerance in soybean. Plant Breed. 2017, 136, 710–719. [CrossRef]27. Chen, W.; Yao, Q.; Patil, G.B.; Agarwal, G.; Deshmukh, R.K.; Lin, L.; Wang, B.; Wang, Y.; Prince, S.J.; Song, L.; et al. Identificationand comparative analysis of differential gene expression in soybean leaf tissue under drought and flooding stress revealed byRNA-Seq. Front. Plant Sci. 2016, 7, 1044. [CrossRef] [PubMed]Agronomy 2021, 11, 92 19 of 2128. Komatsu, S.; Yamamoto, R.; Nanjo, Y.; Mikami, Y.; Yunokawa, H.; Sakata, K. A comprehensive analysis of the soybean genesand proteins expressed under flooding stress using transcriptome and proteome techniques. J. Proteome Res. 2009, 8, 4766–4778.[CrossRef]29. Arora, K.; Panda, K.K.; Mittal, S.; Mallikarjuna, M.G.; Rao, A.R.; Dash, P.K.; Thirunavukkarasu, N. RNAseq revealed the importantgene pathways controlling adaptive mechanisms under waterlogged stress in maize. Sci. Rep. 2017, 7, 10950. [CrossRef]30. Reynoso, M.A.; Kajala, K.; Bajic, M.; West, D.A.; Pauluzzi, G.; Yao, A.I.; Hatch, K.; Zumstein, K.; Woodhouse, M.; Rodriguez-Medina, J.; et al. Evolutionary flexibility in flooding response circuitry in angiosperms. Science 2019, 365, 1291–1295. [CrossRef]31. Staiger, D.; Brown, J.W. Alternative splicing at the intersection of biological timing, development, and stress responses. Plant Cell2013, 25, 3640–3656. [CrossRef]32. Wang, B.B.; Brendel, V. Genomewide comparative analysis of alternative splicing in plants. Proc. Natl. Acad. Sci. USA2006, 103, 7175–7180. [CrossRef] [PubMed]33. Kim, E.; Magen, A.; Ast, G. Different levels of alternative splicing among eukaryotes. Nucleic Acids Res. 2007, 35, 125–131.[CrossRef] [PubMed]34. Kalyna, M.; Simpson, C.G.; Syed, N.H.; Lewandowska, D.; Marquez, Y.; Kusenda, B.; Marshall, J.; Fuller, J.; Cardle, L.; McNicol,J.; et al. Alternative splicing and nonsense-mediated decay modulate expression of important regulatory genes in Arabidopsis.Nucleic Acids Res. 2012, 40, 2454–2469. [CrossRef] [PubMed]35. Drechsel, G.; Kahles, A.; Kesarwani, A.K.; Stauffer, E.; Behr, J.; Drewe, P.; Rätsch, G.; Wachter, A. Nonsense-mediated decay ofalternative precursor mRNA splicing variants is a major determinant of the Arabidopsis steady state transcriptome. Plant Cell2013, 25, 3726–3742. [CrossRef] [PubMed]36. Smith, J.E.; Baker, K.E. Nonsense-mediated RNA decay–A switch and dial for regulating gene expression. BioEssays2015, 37, 612–623. [CrossRef] [PubMed]37. Syed, N.H.; Prince, S.J.; Mutava, R.N.; Patil, G.; Li, S.; Chen, W.; Babu, V.; Joshi, T.; Khan, S.; Nguyen, H.T. Core clock, SUB1,and ABAR genes mediate flooding and drought responses via alternative splicing in soybean. J. Exp. Bot. 2015, 66, 7129–7149.[CrossRef]38. Zou, X.; Jiang, Y.; Zheng, Y.; Zhang, M.; Zhang, Z. Prolyl 4-hydroxylase genes are subjected to alternative splicing in roots ofmaize seedlings under waterlogging. Ann. Bot. 2011, 108, 1323–1335. [CrossRef]39. Van Veen, H.; Vashisht, D.; Akman, M.; Girke, T.; Mustroph, A.; Reinen, E.; Hartman, S.; Kooiker, M.; van Tienderen, P.; Schranz,E.; et al. Transcriptomes of eight Arabidopsis thaliana accessions reveal core conserved, genotype-and organ-specific responses toflooding stress. Plant Physiol. 2016, 172, 668–689. [CrossRef]40. Yasumoto, S.; Terakado, Y.; Matsuzaki, M.; Okada, K. Effects of high water table and short-term flooding on growth, yield, andseed quality of sunflower. Plant Prod. Sci. 2011, 14, 233–248. [CrossRef]41. Van der Merwe, R.; Labuschagne, M.T.; Herselman, L.; Hugo, A. Effect of heat stress on seed yield components and oil compositionin high-and mid-oleic sunflower hybrids. S. Afr. J. Plant Soil 2015, 32, 121–128. [CrossRef]42. Machekposhti, M.F.; Shahnazari, A.; Ahmadi, M.Z.; Aghajani, G.; Ritzema, H. Effect of irrigation with sea water on soil salinityand yield of oleic sunflower. Agric. Water Manag. 2017, 188, 69–78. [CrossRef]43. Hussain, M.; Farooq, S.; Hasan, W.; Ul-Allah, S.; Tanveer, M.; Farooq, M.; Nawaz, A. Drought stress in sunflower: Physiologicaleffects and its management through breeding and agronomic alternatives. Agric. Water Manag. 2018, 201, 152–166. [CrossRef]44. Gao, L.; Lee, J.S.; Hübner, S.; Hulke, B.S.; Qu, Y.; Rieseberg, L.H. Genetic and phenotypic analyses indicate that resistance toflooding stress is uncoupled from performance in cultivated sunflower. N. Phytol. 2019, 223, 1657–1670. [CrossRef] [PubMed]45. Hulke, B.S.; Kleingartner, L.W. Sunflower. In Yield Gains in Major U.S. Field Crops; Smith, S., Diers, B., Specht, J., Carver, B., Eds.;CSSA: Madison, WI, USA, 2014; pp. 433–457. [CrossRef]46. Wample, R.L.; Reid, D.M. Control of adventitious root production and hypocotyl hypertrophy of sunflower (Helianthus annuus) inresponse to flooding. Physiol. Plant 1978, 44, 351–358. [CrossRef]47. Ni, X.L.; Gui, M.Y.; Tan, L.L.; Zhu, Q.; Liu, W.Z.; Li, C.X. Programmed cell death and aerenchyma formation in water-loggedsunflower stems and its promotion by ethylene and ROS. Front. Plant Sci. 2019, 9, 1928. [CrossRef] [PubMed]48. Torres, A.M.; Diedenhofen, U. Sunflower alcohol dehydrogenase genotypes: Germination rates and response to flooding. Environ.Exp. Bot. 1981, 21, 35–44. [CrossRef]49. Jayasekera, G.A.U.; Reid, D.M.; Yeung, E.C. Fates of ethanol produced during flooding of sunflower roots. Can. J. Bot.1990, 68, 2408–2414. [CrossRef]50. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120.[CrossRef]51. Badouin, H.; Gouzy, J.; Grassa, C.J.; Murat, F.; Staton, E.; Cottret, L.; Lelandais-Brière, C.; Owens, G.; Carrere, S.; Mayjonade, B.;et al. The sunflower genome illuminates the evolutionary history of the Asterids and provides new insights into oil metabolismand flowering time. Nature 2017, 546, 148–152. [CrossRef]52. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafastuniversal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [CrossRef]53. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomicfeatures. Bioinformatics 2014, 30, 923–930. [CrossRef]Agronomy 2021, 11, 92 20 of 2154. Rau, A.; Gallopin, M.; Celeux, G.; Jaffrézic, F. Data-based filtering for replicated high-throughput transcriptome sequencingexperiments. Bioinformatics 2013, 29, 2146–2152. [CrossRef]55. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital geneexpression data. Bioinformatics 2010, 26, 139–140. [CrossRef]56. Venables, W.N.; Ripley, B.D. Modern Applied Statistics with S-PLUS; Springer Science & Business Media: New York, NY, USA, 1999.57. Tack, D.C.; Pitchers, W.R.; Adams, K.L. Transcriptome analysis indicates considerable divergence in alternative splicing betweenduplicated genes in Arabidopsis thaliana. Genetics 2014, 198, 1473–1481. [CrossRef]58. Park, E.; Pan, Z.; Zhang, Z.; Lin, L.; Xing, Y. The expanding landscape of alternative splicing variation in human populations. Am.J. Hum. Genet. 2018, 102, 11–26. [CrossRef]59. Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat.Soc. Ser. B Methodol. 1995, 57, 289–300. [CrossRef]60. Bates, D.; Maechler, M.; Bolker, B.; Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 2015, 67, 1–48.[CrossRef]61. Kuznetsova, A.; Brockhoff, P.B.; Christensen, R.H.B. lmerTest package: Tests in linear mixed effects models. J. Stat. Softw. 2017, 82.[CrossRef]62. Fox, J.; Weisberg, S. An R Companion to Applied Regression, 2nd ed.; Sage: Thousand Oaks, CA, USA, 2011; Available online:http://socserv.socsci.mcmaster.ca/jfox/Books/Companion (accessed on 8 May 2019).63. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559.[CrossRef]64. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: Asoftware environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [CrossRef]65. Tian, T.; Liu, Y.; Yan, H.; You, Q.; Yi, X.; Du, Z.; Xu, W.; Su, Z. agriGO v2.0: A GO analysis toolkit for the agricultural community,2017 update. Nucleic Acids Res. 2017, 45, W122–W129. [CrossRef]66. Benjamini, Y.; Yekutieli, D. The control of the false discovery rate in multiple testing under dependency. Ann. Stat.2001, 29, 1165–1188. [CrossRef]67. Supek, F.; Bošnjak, M.; Škunca, N.; Šmuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE2011, 6, e21800. [CrossRef]68. Qi, B.; Yang, Y.; Yin, Y.; Xu, M.; Li, H. De novo sequencing, assembly, and analysis of the Taxodium ‘Zhongshansa’ roots and shootstranscriptome in response to short-term waterlogging. BMC Plant Biol. 2014, 14, 201. [CrossRef]69. Kreuzwieser, J.; Hauberg, J.; Howell, K.A.; Carroll, A.; Rennenberg, H.; Millar, A.H.; Whelan, J. Differential response of graypoplar leaves and roots underpins stress adaptation during hypoxia. Plant Physiol. 2009, 149, 461–473. [CrossRef]70. Sasidharan, R.; Hartman, S.; Liu, Z.; Martopawiro, S.; Sajeev, N.; van Veen, H.; Yeung, E.; Voesenek, L.A.C.J. Signal dynamics andinteractions during flooding stress. Plant Physiol. 2018, 176, 1106–1117. [CrossRef]71. Gasch, P.; Fundinger, M.; Müller, J.T.; Lee, T.; Bailey-Serres, J.; Mustroph, A. Redundant ERF-VII transcription factors bind to anevolutionarily conserved cis-motif to regulate hypoxia-responsive gene expression in Arabidopsis. Plant Cell 2016, 28, 160–180.[CrossRef]72. Xu, X.; Chen, M.; Ji, J.; Xu, Q.; Qi, X.; Chen, X. Comparative RNA-seq based transcriptome profiling of waterlogging response incucumber hypocotyls reveals novel insights into the de novo adventitious root primordia initiation. BMC Plant Biol. 2017, 17, 1–13.[CrossRef]73. Neogy, A.; Garg, T.; Kumar, A.; Dwivedi, A.K.; Singh, H.; Singh, U.; Singh, Z.; Prasad, K.; Jain, M.; Yadav, S.R. Genome-wide transcript profiling reveals an auxin-responsive transcription factor, OsAP2/ERF-40, promoting rice adventitious rootdevelopment. Plant Cell Physiol. 2019, 60, 2343–2355. [CrossRef]74. Grichko, V.P.; Glick, B.R. Ethylene and flooding stress in plants. Plant Physiol. Biochem. 2001, 39, 1–9. [CrossRef]75. Pezeshki, S.R. Wetland plant responses to soil flooding. Environ. Exp. Bot. 2001, 46, 299–312. [CrossRef]76. Glenz, C.; Schlaepfer, R.; Iorgulescu, I.; Kienast, F. Flooding tolerance of Central European tree and shrub species. For. Ecol.Manag. 2006, 235, 1–13. [CrossRef]77. Lee, K.W.; Chen, P.W.; Yu, S.M. Metabolic adaptation to sugar/O2 deficiency for anaerobic germination and seedling growth inrice. Plant Cell Environ. 2014, 37, 2234–2244. [CrossRef]78. Filichkin, S.A.; Priest, H.D.; Givan, S.A.; Shen, R.; Bryant, D.W.; Fox, S.E.; Wong, W.K.; Mockler, T.C. Genome-wide mapping ofalternative splicing in Arabidopsis thaliana. Genome Res. 2010, 20, 45–58. [CrossRef] [PubMed]79. Marquez, Y.; Brown, J.W.; Simpson, C.G.; Barta, A.; Kalyna, M. Transcriptome survey reveals increased complexity of thealternative splicing landscape in Arabidopsis. Genome Res. 2012, 22, 1184–1195. [CrossRef] [PubMed]80. Laloum, T.; Martín, G.; Duque, P. Alternative splicing control of abiotic stress responses. Trends Plant Sci. 2018, 23, 140–150.[CrossRef]81. Reddy, A.S.; Marquez, Y.; Kalyna, M.; Barta, A. Complexity of the alternative splicing landscape in plants. Plant Cell2013, 25, 3657–3683. [CrossRef]82. Soergel, D.A.W.; Lareau, L.F.; Brenner, S.E. Regulation of Gene Expression by Coupling of Alternative Splicing and NMD; Landes Bio-science: Austin, TX, USA, 2013. Available online: https://www.ncbi.nlm.nih.gov/books/NBK6088/ (accessed on 15 May 2019).Agronomy 2021, 11, 92 21 of 2183. Filichkin, S.A.; Mockler, T.C. Unproductive alternative splicing and nonsense mRNAs: A widespread phenomenon among plantcircadian clock genes. Biol. Direct 2012, 7, 20. [CrossRef]84. Matsukura, S.; Mizoi, J.; Yoshida, T.; Todaka, D.; Ito, Y.; Maruyama, K.; Shinozaki, K.; Yamaguchi-Shinozaki, K. Comprehensiveanalysis of rice DREB2-type genes that encode transcription factors involved in the expression of abiotic stress-responsive genes.Mol. Genet. Genom. 2010, 283, 185–196. [CrossRef]85. Peng, X.; Zhang, L.; Zhang, L.; Liu, Z.; Cheng, L.; Yang, Y.; Shen, S.; Chen, S.; Liu, G. The transcriptional factor LcDREB2cooperates with LcSAMDC2 to contribute to salt tolerance in Leymus chinensis. Plant Cell Tissue Organ. Cult. 2013, 113, 245–256.[CrossRef]86. Lee, J.S.; Adams, K.L. Global insights into duplicated gene expression and alternative splicing in polyploid Brassica napus underheat, cold, and drought stress. Plant Genome 2020, e20057. [CrossRef]87. Fukao, T.; Bailey-Serres, J. Plant responses to hypoxia–is survival a balancing act? Trends Plant Sci. 2004, 9, 449–456. [CrossRef][PubMed]88. Bui, L.T.; Novi, G.; Lombardi, L.; Iannuzzi, C.; Rossi, J.; Santaniello, A.; Mensuali, A.; Corbineau, F.; Giuntoli, B.; Perata, P.; et al.Conservation of ethanol fermentation and its regulation in land plants. J. Exp. Bot. 2019, 70, 1815–1827. [CrossRef] [PubMed]89. Loehlin, D.W.; Carroll, S.B. Expression of tandem gene duplicates is often greater than twofold. Proc. Natl. Acad. Sci. USA2016, 113, 5988–5992. [CrossRef]90. Fukao, T.; Xu, K.; Ronald, P.C.; Bailey-Serres, J. A variable cluster of ethylene response factor–like genes regulates metabolic anddevelopmental acclimation responses to submergence in rice. Plant Cell 2006, 18, 2021–2034. [CrossRef]91. Wang, B.; Jin, Q.; Zhang, X.; Mattson, N.S.; Ren, H.; Cao, J.; Wang, Y.; Tao, D.; Xu, Y. Genome-wide transcriptional analysis ofsubmerged lotus reveals cooperative regulation and gene responses. Sci. Rep. 2018, 8, 9187. [CrossRef]92. Yeung, E.; van Veen, H.; Vashisht, D.; Paiva, A.L.S.; Hummel, M.; Rankenberg, T.; Steffens, B.; Steffen-Heins, A.; Sauter, M.;de Vries, M.; et al. A stress recovery signaling network for enhanced flooding tolerance in Arabidopsis thaliana. Proc. Natl. Acad.Sci. USA 2018, 115, E6085–E6094. [CrossRef]93. Razzaque, S.; Elias, S.M.; Haque, T.; Biswas, S.; Jewel, G.N.A.; Rahman, S.; Weng, X.; Ismail, A.M.; Walia, H.; Juenger, T.E.; et al.Gene expression analysis associated with salt stress in a reciprocally crossed rice population. Sci. Rep. 2019, 9, 1–17. [CrossRef]94. Bhattacharjee, A.; Sharma, R.; Jain, M. Over-expression of OsHOX24 confers enhanced susceptibility to abiotic stresses intransgenic rice via modulating stress-responsive gene expression. Front. Plant Sci. 2017, 8, 628. [CrossRef]
Article
eng
Reviewed
Vancouver : University of British Columbia Library
Multidisciplinary Digital Publishing Institute
CC BY 4.0
Faculty
Helianthus annuus L.
crop
abiotic stress
flooding
alcohol dehydrogenase
transcriptomics
alternative splicing
co-expression network
mixed model
fermentation pathway
Genome-Wide Expression and Alternative Splicing in Domesticated Sunflowers (Helianthus annuus L.) under Flooding Stress
Text