The integration of quantitative trait locus mapping and transcriptome studies reveals candidate genes for water stress response in St. Augustinegrass

Abstract

Background

Drought resistance is an increasingly important trait for many plants—including St. Augustinegrass, a major warm-season turfgrass—as more municipalities impose restrictions on frequency and amount of irrigation. Breeding efforts have focused on breeding for drought resistance, and several drought-related quantitative trait loci (QTL) have been identified for St. Augustinegrass in previous studies. However, the molecular basis of this trait remains poorly understood, posing a significant roadblock to the genetic improvement of the species.

Results

This study sought to validate those QTL regions in an independent biparental population developed from two sibling lines, XSA10098 and XSA10127. The drought evaluation in two greenhouse trials showed significant genotype variation for drought stress traits including leaf wilting, percent green cover, relative water content, percent recovery, and the area under the leaf wilting-, percent green cover-, and percent recovery- curves. A linkage map was constructed using 12,269 SNPs, representing the densest St. Augustinegrass linkage map to date. A multiple QTL mapping approach identified 24 QTL including overlapping regions on linkage groups 3, 4, 6, and 9 between this study and previous St. Augustinegrass drought resistance studies. At the transcriptome level, 1965 and 1005 differentially expressed genes were identified in the drought sensitive and tolerant genotypes, respectively. Gene Ontology and KEGG analysis found different mechanisms adopted by the two genotypes in response to drought stress. Integrating QTL and transcriptomics analyses revealed several candidate genes which are involved in processes including cell wall organization, photorespiration, zinc ion transport, regulation of reactive oxygen species, channel activity, and regulation in response to abiotic stress.

Conclusions

By innovatively integrating QTL and transcriptomics, our study advances the understanding of the genetic control of water stress response in St. Augustinegrass, providing a foundation for targeted drought resistance breeding.

Background

Lawn irrigation currently consumes the bulk of home water usage in the US [1]. As more municipalities impose restrictions on irrigation timing and amount, and as climate change leads to more sporadic weather patterns with the potential for more severe drought [2], the development of drought-resistant cultivars is important to ensure turfgrasses will remain a part of landscapes. St. Augustinegrass [Stenotaphrum secundatum (Walt.) Kuntze] is a coarse-textured, warm-season turfgrass commonly grown in the southeastern United States, and it is known for its shade tolerance and weed suppression ability. However, drought resistance in St. Augustinegrass often ranks low among warm-season turfgrasses [3].

In turfgrass species, drought resistance mechanisms are divided into three components: avoidance, escape and tolerance [4]. Drought escape mechanisms involve large-scale changes in growth pattern including varying the timing of the reproductive stage to coincide with more available water or entering dormancy during periods of drought stress and resuming growth when water resources are again available [5]. Because dormancy results in a reduction of turfgrass quality below acceptable levels, it is often a trigger for irrigation among consumers and it may ultimately result in over-watering [6]. Turfgrass breeders are more focused on breeding for genotypes that maintain their green color through avoidance and tolerance. Drought avoidance approaches seek to preventatively address drought through maintenance of plant water content by either increasing the amount of water uptake or decreasing the amount of water loss, thereby avoiding tissue dehydration [5]. Meanwhile, tolerance is the responsive approach of plants to drought stress by sustaining growth despite low water availability and includes osmotic adjustment, maintenance of membrane fluidity, and the accumulation of certain metabolites or proteins [7]. Plants have developed intricate molecular pathways to cope with drought stress, ensuring their survival in water-stress environments. When drought occurs, plants trigger a series of signaling events starting with the perception of water loss, including abscisic acid (ABA) and MAPK [8, 9]. Additionally, drought stress induces the expression of genes involved in osmotic adjustment to help maintain cell turgor [10, 11]. Plants also activate antioxidant systems to combat oxidative stress, which is commonly associated with water stress [12]. Moreover, drought-responsive transcription factors can coordinate a wide range of protective mechanisms, from enhancing root growth to modifying cell wall properties, thereby improving the plant’s ability to retain water and maintain cellular function under drought conditions [13].

Previous studies have established variation in drought resistance among St. Augustinegrass genotypes [6, 14, 15, 16], often noting that polyploids [17] and interploids [15] have better performance. Certain characteristics have been attributed to enhanced drought resistance, namely a deeper root system [16, 18]. However, less work has centered on identifying genomic regions that could be targets for marker-assisted selection (MAS) to improve drought resistance in St. Augustinegrass. Quantitative Trait Loci (QTL) mapping is a method that identifies molecular markers associated with regions controlling traits of interest that can be candidates for MAS. The creation of a high-density linkage map in St. Augustinegrass [19] facilitated the first effort to detect QTL for drought resistance [20]. Using a population developed from the cross between cultivars ‘Raleigh’ and ‘Seville’, 46 significant associations were found, with co-localization of several drought-related traits across greenhouse and field environments [20]. Furthermore, morphological and canopy traits were evaluated in the same population in a later study [21], and several overlapping QTL with those identified in Yu et al. [20] were found. However, as QTL can be specific to a population and influenced by environmental effects, validation of these QTL is necessary before implementing downstream processes such as MAS or genomic selection. In addition, the low-density genetic linkage maps currently available only allow for the identification of QTL regions covering many genes, limiting our ability to identify candidate genes.

High-throughput sequencing technologies give the opportunity to bridge this gap through the integration of QTL mapping and gene expression analysis. Recently, a reference-grade chromosome-scale genome assembly for the popular St. Augustinegrass diploid cultivar ‘Raleigh’ was developed, utilizing data from PacBio CCS, Illumina, and Hi-C technologies [22]. In addition, RNA-Seq has been successfully applied to characterize the transcriptomic profiles in response to abiotic stress [23, 24, 25]. Several studies have shown that integrating QTL mapping and transcriptome study is a robust approach to identifying candidate genes, especially for complex traits [26, 27]. We hypothesize that the integration of QTL mapping and transcriptome analysis will reveal key genetic regions and candidate genes associated with water stress tolerance in St. Augustinegrass. Specifically, we expect to identify novel genes that play a role in physiological processes such as osmotic adjustment, stress-induced signaling, and antioxidant defense, which could serve as targets for improving drought resistance in this species. Thus, the objectives of this study were to: (1) identify and validate drought resistance QTL using a new St. Augustinegrass mapping population; (2) conduct a transcriptomics analysis among two genotypes displaying contrasting performance under drought stress; and (3) through a comprehensive analysis, identify candidate genes within those QTL regions controlling response to drought stress.

Results

Drought-related traits evaluation

The days with the most variance for percent green cover and leaf wilting in 2020 were 17 and 19 days after watering, respectively. In 2021, the days with the most variance for percent green cover and leaf wilting were 14 and 13 days after watering, respectively. The corresponding phenotypic data on those days were selected for further analysis and denoted “_20/21” throughout the rest of the manuscript. For both years, all traits had a wide range of variation and transgressive segregation (Fig. 1). Variance analysis showed that genotype and genotype x year effects were significant for all traits, while year effect was significant for all traits except leaf wilting (Additional file 1). Significant, positive correlations (p-value < 0.001) were found among all traits across years, but correlations were much higher across traits within a year (r = 0.69–0.98 in 2020 and r = 0.74–0.98 in 2021) than between the same trait across years (r = 0.16–0.33) (Fig. 2).

Fig. 1

Histograms of best linear unbiased estimators (BLUEs) for each trait across the two years of evaluation. White triangles indicate the XSA10098 maternal value while black triangles indicate the XSA10127 paternal value. GH = greenhouse, PGC = percent green cover, LW = leaf wilting, PR = percent recovery, RWC = relative water content, AULWC = area under leaf wilting curve, AUPGCC = area under percent green cover curve, AUPRC = area under percent recovery curve

Fig. 2

Pearson correlations for all traits across two years. PGC = percent green cover, LW = leaf wilting, PR = percent recovery, RWC = relative water content, AULWC = area under leaf wilting curve, AUPGCC = area under percent green cover curve, AUPRC = area under percent recovery curve. Color bar indicates correlation coefficient (r). *** indicates p-value < 0.001

SNP discovery and linkage map construction

A total of 940 million raw paired-end reads were obtained from sequencing, and a total of 66,672 SNPs were obtained. Following filtering, 12,269 high-quality SNPs were retained in the final linkage map, which contained nine linkage groups (LGs) corresponding to the nine chromosomes in a St. Augustinegrass haplotype. All linkage groups had between 652 and 2,301 markers with LG 8 having the fewest and LG 9 having the most (Table 1). The length in cM ranged from 77.93 to 173.62 cM. The largest gap between markers was 3.49 cM on LG 5. For most linkage groups, the genomic-informed order produced a shorter linkage group, but for linkage groups 2 and 5, the MDS order proved better and was used in the final map for QTL mapping. Overall, the synteny between markers on the chromosomes and markers on the linkage groups was strong (Additional file 2).

Table 1 Summary of linkage groups including length in cM, total markers, markers/cM, average distance between markers, and the max gap between markers for each linkage group

QTL detection for drought-related traits

Using a stepwise multiple QTL mapping approach, a total of 24 QTL across seven linkage groups were detected (Table 2). Overall, ten QTL were identified across years, 13 QTL were identified from the 2020 trial, and only one QTL was identified from the 2021 trial. For leaf wilting, two QTL were identified across years, three were identified from the 2020 trial, and one was identified from the 2021 trial, with overlapping regions found on LGs 5 and 9. For AULWC, all three QTL, located on LGs 2, 3, 9, were found from the 2020 trial. For PGC, two QTL were identified across years, and four QTL were identified from the 2020 trial, with overlapping regions on LGs 2 and 3. There was one QTL identified for AUPGCC across years and another from the 2020 trial, both of them located at the peak position 2.4 cM on LG 2. For PR, two QTL were identified across years, and one was identified from the 2020 trial, with an overlapping region found on LG 3. Two QTL for AUPRC were identified across years, which were located on LGs 3 and 9. For RWC, one QTL was identified across years, while another was identified from the 2020 trial, with both sharing an overlapping region on LG 3. QTL were also found to be co-localized across different traits in this study. On LG 2, QTL for PGC and AUPGCC were localized on the peak position at 2.4 cM, and a QTL for AULWC was also found in this region. On LG 3, QTL for PR and AUPRC were localized at the peak position 48.31 cM, QTL for AULWC, PGC and RWC were localized at the peak position 101.12 cM, and QTL for LW and RWC were found at the peak position 107.43 cM. In addition, there were a total of 9 QTL for 6 traits overlapping in this interval on LG3. On LG 9, QTL for LW, PR and AUPRC overlapped in the interval between 38.20 cM and 52.05 cM.

Table 2 Significant quantitative trait loci associated with drought-related traits for the XSA10098 X XSA10127 biparental population

QTL overlap across studies

When comparing the results of this study with two previous QTL mapping studies in St. Augustinegrass for drought tolerance traits [20] and morphological characteristics associated with drought tolerance [21], regions of overlap were found on LGs 3, 4, 6, and 9 (Fig. 3). On LG 3, PR-3.1, PR-3.2, PGC-3.1, PGC-3.2, LW-3.1, RWC-3.1, RWC-3.2, and AULWC-3.1 co-localized with RWC-R3.3 from Yu et al. [20]. On LG 4, PGC-4.1 colocalized with GC-R4.1 and Fv/Fm-R4.2 from Yu et al. [20]. On LG 6, LW-6.1 co-localized with LF-R6.1 from Yu et al. [20]. Meanwhile, on LG 9, LW-9.1, LW-9.2, PR-9.1, AUPRC-9.1 and AULWC-9.1 co-localized with qCDr9.1 and qCDs9.1 from Yu et al. [21].

Fig. 3

QTL overlap between the current study and previously identified QTL for drought resistance and morphological characteristics. For black regions from current study: PGC = percent green cover, LW = leaf wilting, PR = percent recovery, RWC = relative water content, AULWC = area under leaf wilting curve, AUPGCC = area under percent green cover curve, AUPRC = area under percent recovery curve. For blue regions from Yu et al. (2019): RWC = relative water content, CC = chlorophyll content, LF = leaf firing, LW = leaf wilting, GC = green cover, NDVI = normalized difference vegetative index. For red regions from Yu et al. (2022): LW = leaf blade width, LL = leaf blade length, CD = canopy density, SGO = shoot growth orientation

Transcriptome assembly

The raw RNA-Seq reads were derived by sequencing cDNA libraries containing leaf and root samples of St. Augustinegrass during water stress. Overall, the RNA sequencing and following quality control process yielded an average of 23,948,230, 24,304,376, 24,560,777, and 24,144,901 clean reads for SWL, SDL, TWL and TDL leaf samples, respectively. For root samples, the quality control process yielded an average of 25,623,156, 24,257,098, 24,462,760, and 23,145,965 clean reads for SWR, SDR, TWR and TDR, respectively. The de novo assembly was used to build St. Augustinegrass transcripts from the clean reads. A total of 997,528 transcripts were obtained with an average length of 651 bp, a maximum length of 55,635 bp, and a minimum length of 201 bp (Table 3). Among all the transcripts, 780,061 unigenes were obtained with an average length of 444 bp. The unigenes were annotated by searching against the public databases (Table 4). The results showed that 291,953 unigenes (37.43%) had significant matches in the Swiss Prot database, 119,598 (15.33%) in the Pfam database, 246,449 (31.59%) in the KEGG database, 289,776 (37.15%) in the GO database, and 75,331 (9.66%) in the EggNOG database. In total, there were 305,220 unigenes (39.13%) successfully annotated in at least one of the public databases, with 26,412 unigenes (3.39%) in all databases.

Table 3 Length of the transcripts and unigenes clustered from the de Novo assembly

Table 4 The annotation of unigenes in different database

Differentially Expressed Genes (DEGs) under water stress

A large number of differentially expressed genes (DEGs) were identified between the different treatment and genotype groups. In this study, we paid particular attention to the DEGs identified for the comparison between with and without drought stress treatment. The comparisons included normal watered vs. drought stress in sensitive genotype (SW vs. SD), and normal watered vs. drought stress in tolerant genotype (TW vs. TD). The number of DEGs identified in the sensitive genotype in response to drought stress was 772 up-regulated and 1193 down-regulated. For the tolerant genotype, 602 DEGs were up-regulated and 403 DEGs were down-regulated in response to drought stress (Fig. 4A, Additional file 3).

Fig. 4

Differentially expressed genes (DEGs) in St. Augustinegrass in response to water stress. (A): Number of DEGs identified in the drought sensitive and tolerant genotypes. (B): Venn diagram of DEGs, showing the unique and overlapping genes expressed among comparison groups. SW vs. SD: normal watered vs. drought stress in sensitive genotype; TW vs. TD: normal watered vs. drought stress in tolerant genotype. Groups A-H indicate different DEGs expression pattern groups

A Venn diagram was used to illustrate the unique and overlapping genes expressed among comparison groups (Fig. 4B). In group A, 118 DEGs were up-regulated in both R and S genotypes. In group B, 180 DEGs were down-regulated in both genotypes. Group C would have contained DEGs that were up-regulated in the S genotype and down-regulated in the T genotype, but none of them fell into this group. Conversely, 2 DEGs that were up-regulated in T but down-regulated in the S genotype fell into group D. Groups E-H were categorized for DEGs that significantly changed only in one genotype and not in the other. Group E included 654 DEGs that were only up-regulated in the S genotype; Group F contained 1011 DEGs that were only down-regulated in S genotype. Group G included 482 DEGs that were only up-regulated in T genotype, and Group H included 223 DEGs that were only down-regulated in T genotype.

GO classification of differentially expressed genes

The DEGs were enriched for GO annotation in terms of molecular function (MF), cellular component (CC), and biological process (BP). In total, 69 MF, 11 CC and 89 BP were classified in group A, while the equivalent number in other groups was: 85/31/123 in group B, 121/34/149 in group E, 154/53/205 in group F, 122/19/177 in group G, and 48/19/90 in group H. Figure 5 shows the top GO classification (gene ratio > 10%) of groups A and B, which include the DEGs commonly up- and down-regulated in both genotypes. Some stress response-related functions, including the abscisic acid metabolic process and ethylene biosynthetic process, are enhanced in both genotypes. In the down-regulated DEGs, several photosynthetic-related functions were shared across genotypes. More details of GO classification of groups A-H can be found in Additional file 4.

Fig. 5

Gene Ontology (GO) classification of the up-regulated DEGs (A) and down-regulated DEGs (B) identified in both genotypes. DEGs were annotated in three categories: biological process, cellular component, and molecular function. Hit% is calculated as the ratio of the input number of DEGs and the number of total annotated genes in this category

KEGG pathway enrichment of differentially expressed genes

The DEGs above were annotated using the KEGG database. In total, 15 pathways were enriched for group A, 33 pathways for group B, 34 pathways for group E, 72 pathways for group F, 15 pathways for group G, and 21 pathways for group H. The top pathway (most number of genes) for group A is 2-oxocarboxylic acid metabolism, the top one for group G is the biosynthesis of secondary metabolites pathway, while the top pathway for groups B, E, F and H is the metabolic pathway. Figure 6 shows the top 20 pathways of groups E-H, which include the DEGs only up- or down-regulated in a single genotype. The MAPK signaling pathway genes were down-regulated in the S genotype, but both up- and down-regulated in the T genotype. The details of the KEGG enrichment are presented in Additional file 5.

Fig. 6

KEGG pathway enrichment of the up-regulated DEGs identified in the S genotype (A) and the T genotype (C) and the down-regulated DEGs identified in the S genotype (B) and the T genotype (D). Hit% is calculated as ratio of input number of DEGs and number of total annotated genes in this pathway

Co-localization of the DEGs and drought QTL

Within the QTL overlapping across studies, twelve DEGs were identified in both the QTL mapping and transcriptomics studies (Table 5; Fig. 7). Among them, two DEGs were identified on chromosome 4, eight DEGs on chromosome 6, and two DEGs on chromosome 9. However, four DEGs encoding probable LRR receptor-like serine/threonine-protein kinase were located in a 16 K bp range on chromosome 6; and they might belong to different exons of the same gene. Two DEGs encoding zinc transporter might also be exons of the same gene. Annotation found these co-localized DEGs represented a wide range of biological processes and pathways including cell wall organization, photorespiration, zinc ion transport, regulation of reactive oxygen species, channel activity, and regulation in response to abiotic stress (Table 5). The DEGs expression patterns show that only TRINITY_DN1864_c1_g1, which encodes a protein cysteine-rich transmembrane module, was up-regulated in the tolerant genotype. Otherwise, all other DEGs were down-regulated in the sensitive genotype, except TRINITY_DN120342_c0_g3, which was down-regulated in both genotypes (Fig. 7).

Table 5 Co-localized differentially expressed gene and overlapping drought QTL regions across the present study

Fig. 7

Expression pattern of differentially expressed genes co-localized in overlapped QTL regions across the present study

Discussion

Drought resistance is a complex trait that is influenced by a number of genes, and typically has low heritability [28]. QTL mapping has been considered as one of the most important approaches for understanding drought genetic architecture [29]. However, low-mapping resolution and population specificity are major constraints in QTL mapping of complex traits. Recent advances in high-throughput genotyping have increased the ability to build high-resolution genetic maps. Compared to the previous linkage maps for St. Augustinegrass developed in Yu et al. [19, 30], the one created in this study represents the densest map to date. More than four times as many markers (12,269) were used in its construction, compared with 2,952 in Yu et al. [19] and 2,257 in Yu et al. [30]. At the same time, the length of all linkage groups ranged from 77.93 to 173.62 cM, differing little from the range of 90.0–197.2 cM in Yu et al. [19] and providing a strong indication that although many more markers were used, the ordering and genotyping was still accurate and did not inflate the cM length.

In addition to the resolution of the genetic map, deciphering and accurately phenotyping drought-related traits has been another challenge for QTL mapping of drought resistance. In our previous studies, we characterized physiological and morphological traits under drought stress, including leaf water content, leaf wilting, leaf firing, NDVI, Fv/Fm, percent green cover, canopy density, leaf length and width, and shoot growth orientation, and successfully mapped QTL for these traits [20, 21]. All these traits were evaluated at the end of drought stress to indicate the final damage on the plants. In this study, we also adopted the area under the curve approach [31], initially developed to measure disease progress, to allow for the evaluation of not only an endpoint in the drought progression but also variation in dry-down rate. Although significant correlations were found for each pair of endpoint and progression traits in this study, and there were overlapping QTL found between them, some specific QTL were also found, indicating progression traits exhibit important value for understanding drought resistance genetics and ultimately improving selection. In addition, the image analysis approach adopted in this study provides stable and accurate evaluation compared to subjective personnel ratings and further supports other studies [32, 33] that the implementation of these approaches increases resolution and efficiency for evaluation of drought tolerance in turfgrass research. Lastly, we also evaluated percent recovery to determine the recovery rate post-stress, which is also an important mechanism that plants use to adapt to drought stress.

Overall, twenty-four QTL were identified from the 2020 and 2021 trials and across years. Surprisingly, despite ample phenotypic variation and a significant genotype effect, only one significant QTL (for leaf wilting) was found in the 2021 trial. In addition to the genetic complexity of drought resistance, environmental variance is known to influence the trait behavior and GxE interaction might also play an important role in its expression, which makes it even more challenging to detect significant QTL and reinforces the need for QTL validation under different conditions. While both experiments were conducted under greenhouse conditions, the maximum temperature was much more stable and higher in 2021, contributing to the faster dry-down. Although there were significant genotype x year interaction effects for all traits, leaf wilting was the only trait that did not show a significant year effect (Additional file 1). That could partially explain why only a LW QTL was detected in the 2021 trial. Likewise, the low correlation between years (r = 0.16– r = 0.33) (Fig. 2) can be adjudicated to inconsistent environmental conditions, which can introduce noise into the data, reducing the power to identify significant QTL associations. While some QTL mapping studies in creeping bentgrass (Agrostis stolonifera L.) have attempted to separate drought [34] and heat conditions [35], when searching for QTL associated with tolerance to those traits, in North Carolina, where our study was conducted, these environmental conditions are often inseparable. Their interaction adds complexity and further demonstrates the need to validate QTL that are stable across environments prior to their use in MAS.

Overlapping QTL emphasized regions on the genome that were responsible for phenotypic variation and well-conserved across the different traits, populations and environments. In this study, overlapping regions across traits and years were found on LG 2, 3 and 9. Of the 24 QTL regions identified in this study, 16 QTL on four LG overlapped with either the previous study on QTL mapping for drought tolerance traits [20] or morphological characteristics associated with drought tolerance [21], which used a different population. Notably, although percent green cover was evaluated in both the current study with 6 identified QTL, and in Yu et al. [20] with 21 identified QTL, only one common region of overlap on LG 4 was identified across the populations. A similar phenomenon occurred for canopy density evaluations in Yu et al. [20, 21]; although the same population was being evaluated, only two common QTL regions were found. These exemplify the difficulty in identifying stable QTL across environments that can be reliably selected for MAS. Thus, the QTL overlapping regions across different populations identified in this study on LGs 3, 4, 6, and 9 will be of interest, which represent strong evidence for their association with the genetic control of St. Augustinegrass’ response to water stress, and as such, the main targets for developing MAS.

The natural outcrossing nature and high-level heterozygosity of St. Augustinegrass limit the fine mapping of large QTL regions. Recently, the integration of QTL mapping and transcriptomics offers a powerful approach to dissecting the genetic basis of complex traits and understanding their underlying molecular mechanisms [27, 36, 37]. By combining genetic mapping with gene expression profiling, researchers can identify candidate genes, unravel regulatory networks, and discover novel genetic variants associated with traits of interest. Among all differentially expressed genes identified in this study, twelve DEGs were found co-localized in QTL overlapping regions (Table 5; Fig. 7). Most of them have been reported to be involved in drought stress response. Probable LRR receptor-like serine/threonine-protein kinases, which is a large gene family that has been implicated for its role in drought stress particularly with abscisic acid signaling [8], were identified. Two zinc transporters were identified on chromosome 6. In a transgenic drought-tolerant maize line, zinc transporter 4 was found to be upregulated compared to the wild-type [38]. An endoglucanase, which is involved in cellulose degradation [39], was annotated on chromosome 9. A transcriptomic analysis found endoglucanases were downregulated in the susceptible line, which was hypothesized to lead to cell wall breakdown and hampered root growth, ultimately contributing to the lack of drought tolerance [40]. Serine hydroxymethyltransferase was down-regulated in the drought sensitive St. Augustinegrass line in this study; similarly, its abundance declined more in a Kentucky bluegrass drought-susceptible line compared to drought-tolerant line [41]. The annotated probable glutathione S-transferase has been widely investigated as an important enzyme that plays a role in ROS scavenging during stress response [12]. Finally, the annotated aquaporin is part of a family of proteins involved in water transport [42]. More interestingly, the DEGs expression pattern shows that only TRINITY_DN1864_c1_g1, which encodes a protein cysteine-rich transmembrane module, was up-regulated in the tolerant genotype. Otherwise, all the other DEGs were down-regulated in S_genotype, except TRINITY_DN120342_c0_g3 was down-regulated in both genotypes (Fig. 7), indicating that the expression change of these genes might partially contribute to less drought resistance in the sensitive genotype.

Beyond the DEG co-localized in QTL regions, the transcriptomics study provided an abundance of gene expression profiles and metabolic pathways of drought tolerant and sensitive genotypes, which could help us to understand how their genetic differences manifest at the level of gene expression, and how these differences contribute to differential performance under drought stress. The MAPK cascade is considered a major signal transducer, playing a vital role in drought stress, generally by responding to ABA and regulating ROS production [9]. Numerous components of MAPK cascades have been reported to respond to drought in crops. Recent RNA-Seq studies found that the expression of several MAPK transcripts changed under drought stress in rice, wheat, cotton and maize, highlighting the importance of MAPKs in drought [43, 44, 45, 46]. In St. Augustinegrass, we found that the MAPK signaling pathway was enriched in down-regulated genes for both genotypes, but it was only enriched in up-regulated genes for the tolerant genotype (Fig. 6, Additional file 5). The up-regulation of MAPK signaling pathway genes -- including five PP2C genes (Probable protein phosphatase 2 C)-- were reported to function as a switch at the center of the ABA signaling network under osmotic stress conditions [47], and might have contributed to the superior performance of the tolerant genotype during drought. In addition, we also noticed another pathway, starch and sucrose metabolism, which was enriched in up-regulated genes in both genotypes but only enriched in down-regulated genes for the sensitive genotype (Fig. 6, Additional file 5). Starch and sucrose are key molecules in mediating plant responses to abiotic stresses. Plants generally remobilize starch to provide energy and carbon at times when photosynthesis may be potentially limited during drought stress. The released sugars and other derived metabolites support plant growth under stress and function as osmoprotectants and compatible solutes to mitigate the negative effects of the stress [10, 11]. In the sensitive genotype, downregulation of key genes -- including Beta-glucosidase 31 (TRINITY_DN18016_c0_g1) and Sucrose synthase 4 (TRINITY_DN27451_c0_g1)– was observed. These genes have been previously implicated in drought stress response via starch and sucrose metabolism pathways [48, 49]. We speculate that the normal starch and sucrose metabolism was more affected by drought stress in sensitive genotype, which might lead to its poor performance under stress. Finally, when we investigated the DEG showing the opposite changes in the two genotypes, we found two genes TRINITY_DN44494_c0_g1 (Protein DETOXIFICATION 19) and TRINITY_DN4181_c0_g1 (Probable cytokinin riboside 5’-monophosphate phosphoribohydrolase LOGL9) were up-regulated in the tolerant genotype but down-regulated in the sensitive genotype (Additional file 3). Protein DETOXIFICATION has been reported to be of significance in the translocation of abscisic acid (ABA), a phytohormone with profound role in plants under various abiotic stress conditions. Arabidopsis lines over-expressing a cotton DETOXIFICATION gene were highly tolerant to drought, salt, and cold stress with high production of antioxidant enzymes and significantly reduced levels of oxidants [50]. The LONELY GUY (LOG) gene was reported to be involved in biosynthesis of cytokinins, which are generally considered to be negative regulators of stress [51]. In Arabidopsis, the expression levels of GhLOG were changed by drought and salt stresses, and the overexpression of GhLOG3 improved salt tolerance potentially though regulation of root growth [52]. These pathways and genes identified from transcriptomics analysis might partially explain the different drought tolerance levels of the two genotypes, pointing to their potential value for marker assisted selection of St. Augustinegrass. However, the function of these genes involved in drought response still needs to be further investigated and validated.

Conclusion

Although drought resistance is a complex mechanism in plants, the integration of different genetic methods may help improve our understanding of how plants adapt to drought stress. Overall, our study identified and validated drought resistance QTL on different populations in beneficial from higher resolution linkage map and higher throughput phenotyping approaches. In addition, the comparative transcriptomic analysis between two contrasting lines identified differentially expressed genes and pathways. Finally, our study integrated St. Augustinegrass transcriptomic resources together with QTL mapping to explore candidate genes in response to drought stress. Our work has implications both for the scientific understanding of drought stress in grasses as well as for practical applications aimed at breeding St. Augustinegrass with superior drought resistance.

Materials and methods

Mapping population development

For QTL mapping and validation, a pseudo-F₂ mapping population was derived from the cross of North Carolina State University St. Augustinegrass breeding lines XSA10098 and XSA10127, which are both full-sib progeny of crosses between cultivar ‘Raleigh’ and ‘Seville’. These breeding lines were selected based on their contrasting responses to drought under both greenhouse and field conditions from Yu et al. [20]. Crosses were verified based on the segregation at two SNP loci using allele-specific primers and the PCR Allele Competitive Extension (PACE) assay (3CR Bioscience, Harlow, United Kingdom).

Greenhouse drought evaluations

The mapping population was vegetatively propagated into 15 cm (diameter) x 11 cm (deep) pots containing a 1:3 mix of sand and Fafard potting mix (Conrad Fafard Inc., Agawam, MA). Plants were allowed to grow for 12 weeks under normal watering conditions to fully establish. The progenies, along with the parents, were organized into a randomized complete block design with three replicates. The first trial (GH20) took place from 3 December 2020 to 21 December 2020 and lasted 19 days. The second trial (GH21) took place from 26 April to 9 May and lasted for 14 days. In the beginning of each trial, pots were watered to capacity, and water was subsequently withheld to initiate drought. The dry-down concluded when each replicate had at least 60% of plants showing leaf wilting symptoms, then watering was initiated once again for 10 days of recovery.

Leaf wilting (LW) ratings were taken daily on a 1 to 5 scale with 5 indicating a healthy, symptom-free turf stand, while a score of 1 indicated desiccated turf that had lost all its color (Fig. 8). On the last day of the dry-down, leaf tissue was collected to determine leaf relative water content (RWC) using the following equation,

Fig. 8

Leaf wilting rating scale for St. Augustinegrass experiencing drought stress. A score of 5 indicates healthy, symptom-free turf while a score of 1 indicates desiccated turf that has lost all its green color

$$\:RWC=\frac{FW-DW}{TW-DW}\:x\:100$$

where FW is fresh weight, DW is dry weight, and TW is turgid weight as described in Yu et al. [20]. A panoramic 180i game camera (Moultrie Feeders, LLC, Birmingham, AL) was mounted above each replicate and set to collect an image of plants prior to drought stress, and once a day throughout dry-down and recovery. Percent green cover was calculated from these images in ArcMap (ESRI, Redlands, CA). Briefly, images were loaded, and a training set was developed to classify green and brown pixels. A supervised classification method analyzed the images and individual percent green cover values were extracted for each pot and normalized to the first date to account for any differences in initial establishment. The values in dry-down were defined as the trait percent green cover (PGC), and values in recovery period were defined as the trait percent recovery (PR). In addition, the area under the progress curve for leaf wilting, percent green cover, and percent recovery (referred to as AULWC, AUPGCC, and AUPRC, respectively) was calculated using the agricolae R package [53].

Statistical analysis

A fully random model using both entry and replicate as random effects was fit for each day of percent green cover and leaf wilting to calculate variances. The day with the most entry variance for percent green cover and leaf wilting from each year was selected to be used for further analysis. Pearson correlations between all traits were generated using the corrplot package in R [54].

For variance analysis, a mixed model was first fit for each trait across years with genotype as fixed and year and replicate as random terms (Eq. 1.) using ASReml-R package [55]

$$\:y\:=\:\mu\:1\:+{{X}_{1}g+Z}_{1}r\:+{Z}_{2}y+{Z}_{3}gy+e$$

(Eq. 1)

where y is the vector of phenotypic values, $\:\mu\:$ is the overall mean, 1 is a vector of ones, X and Z represent the incidence matrices for fixed and random effects, respectively; g is the fixed vector of genotype effects, and r, y and gy are the random vector of effects of rep, year and the interaction between genotype and year, respectively. Using a Wald test for fixed effects and a likelihood ratio test (LRT) for random effects, the significance of each term was tested, and variance components were estimated with restricted maximum likelihood. Best linear unbiased estimators (BLUEs) generated from Eq. 1 were used for QTL analysis of traits across years.

For QTL analysis within each year, a mixed model was fit for each trait within a year with entry as a fixed term and rep as random (Eq. 2.)

$$\:y\:=\:\mu\:1\:+{{X}_{1}g+Z}_{1}r+e$$

(Eq. 2)

with terms the same as in Eq. 1. Best linear unbiased estimators were generated for each entry using Eq. 2.

Library construction and sequencing

Young leaves were collected from each of the progeny as well as the parents (five replicates per parent to improve coverage), and genomic DNA was extracted using the CTAB method modified from Afanador et al. [56]. DNA quality and purity was validated using a 1% agarose gel, and DNA concentration was quantified using a Hoefer DQ 300 fluorometer (Hoefer Scientific Instruments, San Francisco, CA). The GBS library was prepared by the University of Wisconsin Biotechnology Center DNA Sequencing Facility by adapting the methods in Elshire et al. [57]. Briefly, PstI and MspI (New England Biolabs, Ipswich, MA) were used to digest 100 ng of DNA, and barcoded adapters were ligated using T4 ligase (New England Biolabs, Ipswich, MA). Samples were pooled and amplified, and SPRI beads were used to remove adapter dimers. The library was quantified using Qubit^® dsDNA HS Assay Kit (Life Technologies, Grand Island, NY), and quality was confirmed using the Agilent Tapestation (Agilent Technologies Inc., Santa Clara, CA). Sequencing was done on an Illumina NovaSeq6000 (Illumina, San Diego, CA), at the University of Wisconsin Biotechnology Center DNA Sequencing Facility.

SNP calling

The v.4.1 GBS-SNP-CROP (GSC) pipeline [58] was used in conjunction with FASTQ [59] to process raw sequencing data and call SNPs. The GSC pipeline identified reads with barcode sequences next to a cut site, trimming them and retaining only those reads. Fastp was used to trim any remaining barcodes and poly-G sequences and reads that had a length of at least 32 bases and at least 40% of bases with a Phred quality score of 30 or greater were retained. Sequences were demultiplexed and aligned to the St. Augustinegrass reference genome [22], using the GSC pipeline. Five samples of each parental genotype were combined. Filtering and genotypic calls were made using the following parameters: mnHoDepth0 = 3, mnHoDepth1 = 10, mnHetDepth = 2, altStrength = 0.96, mnAlleleRatio = 0.05, mnAvgDepth = 3, and mxAvgDepth = 300. SNPs with less than 70% genotype calls were removed.

Linkage map construction

The R package MAPpoly (v. 0.3.0) was used to jointly construct the linkage map [60] considering markers of all possible segregation types. Briefly, a chi-square test using a Bonferroni correction to achieve a global significance level of 5% was used to assess segregation distortion, and any markers deviating from the expected ratios (1:2:1 or 1:1) were excluded. The pairwise recombination fractions were calculated between all markers and considering all possible linkage phases. The resulting recombination fraction matrix based on the most likely linkage phases was used for grouping and ordering markers inside linkage groups. The unweighted pair-group method with arithmetic mean (UPGMA) clustering algorithm was used to form linkage groups, and the multidimensional scaling (MDS) [61] algorithm was used to order markers. The reference genome was used as supporting information to group markers; thus linkage groups were formed only by markers that agreed with both the UPGMA and reference genome information. A hidden Markov model was used to assess the best linkage phases and re-estimate the recombination fractions between markers for all linkage groups, considering the orders provided by the MDS algorithm and the reference genome. The best order for each linkage group was selected based on the pattern of the recombination frequency heatmap, the resulting map size, and the number of retained markers. A global error rate of 5% was considered to account for sequencing and genotyping errors.

QTL mapping

The R package QTLpoly (v. 0.2.1), which was designed for outcrossing species, was used to conduct QTL mapping [62], along with BLUEs for each trait from each year and from the joint analysis across years. Conditional genotype probabilities were first calculated at a step of 1 centimorgan (cM) using the “calc_genoprob” function in MAPpoly. The random-effect multiple interval mapping (REMIM) model was used to build a multiple QTL model for each trait [62]. First, a null model with no QTL was fit, then QTL were added in a forward step one by one until no more QTL reached the threshold (P < 0.01). A backward elimination was performed with a more stringent threshold level (P < 0.0001) to remove QTL that presented non-significant effects conditional on the other ones. The process was repeated until no QTL was added nor dropped from the model, then a last round of model refinement was performed to adjust QTL position and effect estimates conditional on the remaining QTL in the final model.

QTL validation with previous studies

The QTL in this study were compared to those identified in the St. Augustinegrass biparental mapping population ‘Raleigh’ x ‘Seville’ in Yu et al. [20], which analyzed drought tolerance traits, and those in Yu et al. [21], which analyzed morphological characteristics related to drought tolerance. The comparison was made possible by mapping the SNP marker sequences on ‘Raleigh’ x ‘Seville’ maps to the reference genome to obtain physical positions. Regions of chromosomes that contained co-localized QTL across studies were visualized in MapChart [63].

Materials and treatment for transcriptomics study

XSA10098 (tolerant) and ‘Raleigh’ (sensitive) showing contrasting drought tolerance levels were used in the transcriptomics study [20]. Pots were established and maintained as the above mapping population in the greenhouse, and then subjected to water stress by cutting off irrigation. Fully expanded leaves with similar growing stages were collected when the drought sensitive line showed leaf wilting symptoms (14 days after cutting off irrigation). Three biological replicates were included for each genotype. Samples were frozen in liquid nitrogen immediately and stored at -80 ^oC until RNA extraction. The abbreviations [S/T][W/D][L/R] [1, 2, 3] were used for sample names in this study: the first letter indicates genotype (sensitive or tolerant), the second letter indicates treatment (watered or drought), the third letter indicates tissue (leaf or root), and the number indicates replication.

RNA extraction, library Preparation and sequencing

RNA extractions and the sequencing library were prepared according to Brown et al. [25]. Briefly, the total RNA was extracted using RNeasy Plant Mini Kit (Qiagen, Germany), and then subjected to quality and quantity checks (RIN value ≥ 7; total amount ≥ 1.5 ug). For library preparation, mRNA was enriched via oligo(dT) beads and fragmented in fragmentation buffer. First strand cDNA was synthesized using random hexamers and reverse transcriptase, then second strand synthesis buffer (Illumina, San Diego, CA) was added along with dNTPs, RNase H and Escherichia coli polymerase I for second strand synthesis. After a round of purification, terminal repair and poly-A tailing, Illumina PE adapters were ligated and 200 bp cDNA fragments were preferentially selected. Then, cDNA fragments with ligated adapters were enriched using PCR primers. The library quality was assessed by quantifying with a Qubit 2.0 fluorometer and qPCR, and insert size was checked on an Agilent 2100. The library was sequenced using an Illumina NovaSeq sequencing platform and 150 bp paired end reads were generated.

Quality control, transcriptome assembly and gene annotation

The raw reads were processed to remove reads containing adapters, uncertain nucleotides (N > 10%), and low-quality reads. Obtained clean reads were evaluated on Q20, Q30, GC content, and used for downstream analyses. Transcriptome assembly was built from paired-end reads from leaf samples in this study and root samples generated in a separate study [64] using Trinity with min_kmer_cov set to 2 by default and with all other parameters kept as default [65]. Gene function was annotated against the Swiss-Prot and Pfam databases and assigned to functional categories in the EggNOG, GO and KEGG databases by searching BLASTx and BLASTp with an E value cutoff 10^− 5.

Differential expression analysis

To explore the gene expression change during water stress in leaf samples, the read counts were adjusted by edgeR package for gene quantification [66]. Differential gene expression analysis was performed between samples with and without water stress for each genotype. Differential gene expression analysis of the samples was performed using DESeq R package [67], and the significant differential expression was filtered using p value ≤ 0.001 and| log₂(Fold Change)| ≥ 1.

GO and KEGG enrichment analysis

Differentially expressed genes (DEGs) were assigned to Gene Ontology (GO) and KEGG pathway enrichment analysis via the GOseq R package [68]. Gene length bias was adjusted by calculating a probability weighting function (PWF). A P value < 0.05 was used as the significance threshold.

Co-Localization of the DEGs and drought QTL

Sequences of DEGs were aligned to the reference genome to obtain their physical positions. Then, the respective genomic regions of the above co-localized QTL were compared to the locations of the identified DEGs.