Integration of QTL and Transcriptome Studies Reveals Candidate Genes for Water Stress Response in St. Augustinegrass

doi:10.21203/rs.3.rs-4701704/v1

Download PDF

Research Article

Integration of QTL and Transcriptome Studies Reveals Candidate Genes for Water Stress Response in St. Augustinegrass

https://doi.org/10.21203/rs.3.rs-4701704/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

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 QTL have been identified for St. Augustinegrass in our previous studies. However, the molecular basis of this trait is still less understood, which has been a significant roadblock for 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 from a total of 12,269 SNPs, representing the densest St. Augustinegrass linkage map to date. Twenty-four QTL were identified from a multiple QTL mapping approach, and overlapping regions from this study and previous St. Augustinegrass drought resistance studies were found on linkage groups 3, 4, 6, and 9. 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

These results represent a step toward understanding the genetic control of water stress response in St. Augustinegrass and provide a theoretical basis for genetic improvement of drought resistance in this species.

QTL mapping

transcriptome

drought stress

St. Augustinegrass

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].

Previous studies have established variation in drought resistance among St. Augustinegrass genotypes [6, 8–10], often noting that polyploids [11]and interploids [9] have better performance. Certain characteristics have been attributed to enhanced drought resistance, namely a deeper root system [10, 12]. 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 [13] facilitated the first effort to detect QTL for drought resistance [14]. 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 [14]. Furthermore, morphological and canopy traits were evaluated in the same population in a later study [15], and several overlapping QTL with those identified in [14] 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 [16]. In addition, RNA-Seq has been successfully applied to characterize the transcriptomic profiles in response to abiotic stress [17–19]. Several studies have shown that integrating QTL mapping and transcriptome study is a robust approach to identifying candidate genes, especially for complex traits [20, 21]. The differentially expressed genes from transcriptomics studies can be used in combination with identified QTL regions to narrow down the number of candidate genes. Thus, the objectives of this study were be to 1) develop a new St. Augustinegrass mapping population segregating for drought resistance and to validate the previously identified drought resistance QTL; 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.

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

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–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.**
Linkage Group	Length (cM)	Number of markers	Markers/cM	Average distance between markers (cM)	Max Gap (cM)
1	101.89	1275	12.51	0.08	2.2
2	97.73	1313	13.43	0.07	2.05
3	117.39	1444	12.30	0.08	2.76
4	102.34	1579	15.43	0.06	2.23
5	129.28	1331	10.30	0.10	3.49
6	111.32	1119	10.05	0.10	2.12
7	94.79	1255	13.24	0.08	1.88
8	77.93	652	8.37	0.12	3.30
9	173.62	2301	13.25	0.08	3.12

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 to 52.05 cM.

Table 2

**Significant quantitative trait loci associated with drought-related traits for the XSA10098 x XSA10127 biparental population.**
QTL	Trait	LG	Peak Position (cM)	% Var	Confidence Interval	Peak Physical Position (bp)	Physical Position Range (bp)
LW-5.1	LW	5	23.66	15.55	16.68–36.06	3363900	2272746 − 4598389
LW-9.2	LW	9	45.27	25.61	37.37–59.13	9573540	7294000–11467984
LW-3.1	LW_20	3	107.43	19.15	101.12–108.19	45180406	43879062–45318615
LW-6.1	LW_20	6	32.35	9.6	20.57–52.29	10814818	7656452–28621279
LW-9.1	LW_20	9	43.19	22.5	38.20–59.13	9326977	7892898–11467984
LW-5.2	LW_21	5	30.31	17.04	16.68–44.13	4161518	2272746 − 8515041
AULWC-2.1	AULWC_20	2	0	28.26	0.00–5.01	8589918	8589918–10251920
AULWC-3.1	AULWC_20	3	101.12	14.92	98.24–108.19	43879062	43307483–45318615
AULWC-9.1	AULWC_20	9	7.11	11.53	7.11–51.01	2791778	2791593–10550638
PGC-2.2	PGC	2	10.24	24.53	0–26.3	11808394	8589918–17037003
PGC-3.2	PGC	3	101.12	12.65	48.31–108.19	43879062	21153247–45318615
PGC-1.1	PGC_20	1	5.31	13.11	5.31–8.71	13405134	13405134–14189176
PGC-2.1	PGC_20	2	2.4	37.06	2.40–3.08	9330187	9330187–9524355
PGC-3.1	PGC_20	3	99.23	12.02	87.32–107.43	43501307	40644353–45180406
PGC-4.1	PGC_20	4	42.27	7.82	34.22–71.06	14146107	9601868–33495704
AUPGCC-2.1	AUPGCC	2	2.4	26.29	0–10.24	9330187	8589918–11808394
AUPGCC-2.2	AUPGCC_20	2	2.4	30.69	0.00–6.12	9330187	8589918–10906196
PR-3.1	PR	3	48.31	14.26	37.06–98.24	21153247	18949946–43307483
PR-9.1	PR	9	43.19	29.08	38.20–54.21	9326977	7892898–11027839
PR-3.2	PR_20	3	98.24	18.3	48.31–108.19	43307483	21153247–45318615
AUPRC-3.1	AUPRC	3	48.31	16.28	41.08–52.54	21153247	19612007–21898037
AUPRC-9.1	AUPRC	9	43.19	27.29	37.37–52.05	9326977	7294000–10721565
RWC-3.1	RWC	3	101.12	15.93	48.31–109.59	43879062	21153247–45704345
RWC-3.2	RWC_20	3	107.43	23.31	96.68–108.19	45180406	42977450–45318615
LG = linkage group, 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. Each trait was analyzed using across years, 2020 and 2021 trial.

QTL overlap across studies

When comparing the results of this study with two previous QTL mapping studies in St. Augustinegrass for drought tolerance traits [14] and morphological characteristics associated with drought tolerance [15], 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. [14]. On LG 4, PGC-4.1 colocalized with GC-R4.1 and Fv/Fm-R4.2 from Yu et al. [14]. On LG 6, LW-6.1 co-localized with LF-R6.1 from Yu et al. [14]. 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. [15].

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 23948230, 24304376, 24560777, and 24144901 clean reads for SWL, SDL, TWL and TDL leaf samples, respectively. For root samples, the quality control process yielded an average of 25623156, 24257098, 24462760, and 23145965 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 997528 transcripts were obtained with an average length of 651 bp, a maximum length of 55635 bp, and a minimum length of 201 bp (Table 3). Among all the transcripts, 780061 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 291953 unigenes (37.43%) had significant matches in the Swiss Prot database, 119598 (15.33%) in the Pfam database, 246449 (31.59%) in the KEGG database, 289776 (37.15%) in the GO database, and 75331 (9.66%) in the EggNOG database. In total, there were 305220 unigenes (39.13%) successfully annotated in at least one of the public databases, with 26412 unigenes (3.39%) in all databases.

Table 3

**Length of the transcripts and unigenes clustered from the de novo assembly**
Category	Transcript	Unigene
Total	997528	780061
Max length	55635 bp	55635 bp
Min length	201 bp	201 bp
Average length	651 bp	444 bp
200–500 bp	697409	633059
500–1000 bp	136590	87389
1000–2000 bp	97721	43018
> 2000 bp	65808	16595
N50	1207 bp	590 bp
N90	251 bp	128 bp
Note: N50 is defined as the length of the shortest contig at 50% of the total assembly length; N90 was counted in the similar way.

Table 4

**The annotation of unigenes in different database**
	No. Unigenes	Percentage (%)
Annotated in Swiss Prot	291953	37.43
Annotated in Pfam	119598	15.33
Annotated in KEGG	246449	31.59
Annotated in GO	289776	37.15
Annotated in EggNOG	75331	9.66
Annotated in all databases	26412	3.39
Annotated in at least one database	305220	39.13
Total unigenes	780061

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

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.

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.

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 16K 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.**
Gene_id	Protein	Chr.	Physical Position (bp)	Function annotation
TRINITY_DN10664_c0_g1	Aquaporin NIP2-2	4	10146118–10146422	Channel activity, transmembrane transporter activity
TRINITY_DN119686_c0_g1	Probable glutathione S-transferase DHAR2, chloroplastic	4	10432628–10432944	Involved in scavenging reactive oxygen species (ROS)
TRINITY_DN7546_c0_g2	Probable LRR receptor-like serine/threonine-protein kinase	6	7942815–7942890	Involved in abiotic stress responses; induction of the activities of antioxidative enzymes; control of stomatal development
TRINITY_DN7546_c2_g1		6	7944167–7944464
TRINITY_DN17100_c0_g2		6	7957142–7957310
TRINITY_DN85968_c0_g1		6	7957612–7958072
TRINITY_DN1864_c1_g1	Protein cysteine-rich transmembrane module	6	8148983–8149282	Regulation of response to salt stress; regulation of reactive oxygen species biosynthetic process
TRINITY_DN283583_c0_g1	Retrovirus-related Pol polyprotein from transposon RE1	6	8247717–8247988	DNA integration; DNA recombination
TRINITY_DN8017_c1_g1	Zinc transporter	6	8283935–8284983	Zinc ion transmembrane transport
TRINITY_DN213079_c0_g3	Zinc transporter	6	8285892–8286045	Zinc ion transmembrane transport
TRINITY_DN12130_c0_g1	Serine hydroxymethyltransferase, mitochondrial	9	5977060–5977249	Involved in photosynthetic pathway, in response to drought and salt stress
TRINITY_DN120342_c0_g3	Endoglucanase 10	9	6093599–6093797	Cell wall organization
Chr.: chromosome.

Drought resistance is a complex trait that is influenced by a number of genes, and typically has low heritability [22]. QTL mapping has been considered as one of the most important approaches for understanding drought genetic architecture [23]. 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. [13, 24], 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. [13] and 2,257 in Yu et al. [24]. At the same time, the length of all linkage groups ranged from 77.93–173.62 cM, differing little from the range of 90.0–197.2 cM in Yu et al. [13] 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 [14, 15]. 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 [25], 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 [26, 27]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 [28] and heat conditions [29], 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 [14] or morphological characteristics associated with drought tolerance [15], 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. [14] 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. [14, 15]; 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 [21, 30, 31]. 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 [31], 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 [32]. An endoglucanase, which is involved in cellulose degradation [33], 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 [34]. 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 [35]. The annotated probable glutathione S-transferase has been widely investigated as an important enzyme that plays a role in ROS scavenging during stress response [36]. Finally, the annotated aquaporin is part of a family of proteins involved in water transport [37]. 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 [38]. 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 [39–42]. 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-regulated MAPK signaling pathway genes might have contributed to better 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 [43, 44]. 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 [45]. The LONELY GUY (LOG) gene was reported to be involved in biosynthesis of cytokinins, which are generally considered to be negative regulators of stress [46]. 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 [47]. 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.

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 integrated newly developed St. Augustinegrass transcriptomic resources together with QTL mapping on different populations to validate previous identified QTL and 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.

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. [14]. 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,

$$\: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. [14]. 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 [48].

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 [49].

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 [50]

$\: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. [51]. 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. [52]. 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 [53]was used in conjunction with fastp [54] 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 [16], 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 [55] 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) [56] 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 [57], 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 [57]. 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.

Materials and treatment for transcriptomics study

XSA10098 (tolerant) and ‘Raleigh’ (sensitive) showing contrasting drought tolerance levels were used in the transcriptomics study [14]. 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–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. [19]. 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 [58] using Trinity with min_kmer_cov set to 2 by default and with all other parameters kept as default [59]. 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 [60]. 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 [61], 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 [62]. Gene length bias was adjusted by calculating a probability weighting function (PWF). A P value < 0.05 was used as the significance threshold.

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. [14], which analyzed drought tolerance traits, and those in Yu et al. [15], 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].

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.

AULWC

area under the leaf wilting curve

AUPGCC

area under the percent green cover curve

AUPRC

area under the percent recovery curve

BLAST

basic local alignment search tool

BLUE

best linear unbiased estimator

DEG

differentially expressed gene

gene ontology

GSC

GBS-SNP-CROP

linkage group

LRT

likelihood ratio test

leaf wilting

MAS

marker-assisted selection

MDS

multidimensional scaling

MQM

multiple QTL mapping

PGC

percent green cover

percent recovery

QTL

quantitative trait locus/loci

RWC

relative water content

SNP

single nucleotide polymorphism

UPGMA

unweighted pair-group method with arithmetic mean.

Availability of data and materials

All the raw sequencing data were deposited into the database of sequence read archive (SRA) in National Center for Biotechnology Information (NCBI) under BioProject ID: PRJNA1134789. The other datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

GR: Investigation, Formal Analysis, Visualization, Writing - Original Draft & Review

XY: Conceptualization, Project Administration, Investigation, Formal Analysis, Visualization, Writing - Original Draft & Review & Editing

GG: Methodology, Software, Methodology, Software, Writing - Review

SG: Validation, Formal Analysis

AD: Formal Analysis

BG: Data Curation, Formal Analysis

AS: Methodology

AH: Methodology

SM: Conceptualization, Resources, Supervision, Funding Acquisition, Writing - Review & Editing.

All the authors reviewed the manuscript.

Funding

This work was supported by a Specialty Crop Research Initiative grant [2019-51181-30472] from the USDA National Institute for Food and Agriculture as well as funding from the North Carolina State University Center for Turfgrass Environmental Research and Education.

Acknowledgments

We thank the University of Wisconsin Biotechnology Center DNA Sequencing Facility (Research Resource Identifier – RRID:SCR_017759) for providing GBS library preparation and sequencing services.

Environmental Protection Agency. Water efficiency management guide landscaping and irrigation. EPA. 2017. 832-F-17-016b. https://www.epa.gov/sites/default/files/2017-12/documents/ws-commercialbuildings-waterscore-irrigation-landscape-guide.pdf.
Kunkel KE, Stevens LE, Stevens SE, Sun L, Janssen E, Wuebbles D, et al. Regional climate trends and scenarios for the U.S. national climate assessment. Part 2. Climate of the Southeast U.S. NOAA Technical Report NESDIS.142-2. 2013; https://nesdis-prod.s3.amazonaws.com/migrated/NOAA_NESDIS_Tech_Report_142-2-Climate_of_the_Southeast_US.pdf.
Harivandi MA, Baird J, Hartin J, Henry M, Shaw D. (2009). Managing turfgrasses during drought. University of California Division of Agriculture and Natural Resources. 2009; Publication 8395. https://doi.org/10.3733/ucanr.8395.
Levitt J. Responses of plants to environmental stress. 2nd ed. Academic Press, New York; 1980.
Huang B. Mechanisms and strategies for improving drought resistance in turfgrass. Acta Horticulturae, 2008;783:221-7.
Steinke K, Chalmers D, Thomas J, White R, Fipps G. Drought response and recovery characteristics of St. Augustinegrass cultivars. Crop Science, 2010;50:2076-83.
Nilsen ET, Orcutt DM. The physiology of plants under stress. Wiley, New York; 1996.
Katuwal KB, Jespersen D, Bhattarai U, Chandra A, Kenworthy KE, Milla-Lewis SR et al. Multi-locational screening identifies new drought tolerant warm-season turfgrasses. Crop Science, 2022;62:1614-30.
Meeks M, Chandra A. Drought response and minimal water requirements of diploid and interploid St. Augustinegrass under progressive drought stress. Crop Science, 2020;60:1048-63.
Zhang J, Poudel B, Kenworthy K, Unruh JB, Rowland D, Erickson JE, Kruse J. Drought responses of above-ground and below-ground characteristics in warm-season turfgrass. Journal of Agronomy and Crop Science, 2019;205:1-12.
Busey P. St. Augustinegrass, Stenotaphrum secundatum (Walt.) Kuntze. In: Casler MD, Duncan RR, editors. Turfgrass biology, genetics, and breeding. John Wiley & Sons, Hoboken, NJ; 2003.
Chandra A, Genovesi AD, Wherley BG, Metz SP, Reinert JA, Wu YZ et al. Registration of ‘DALSA 0605’ St. Augustinegrass. Journal of Plant Registrations, 2015;9:27-34.
Yu X, Kimball JA, Milla-Lewis SR. High density genetic maps of St. Augustinegrass and applications to comparative genomic analysis and QTL mapping for turf quality traits. BMC Plant Biology, 2018;18:1-13.
Yu X, Brown JM, Graham SE, Carbajal EM, Zuleta MC, Milla-Lewis SR. Detection of quantitative trait loci associated with drought tolerance in St. Augustinegrass. PLoS ONE, 2019;14: 1-14.
Yu X, Lara NAH, Carbajal EM, Milla-Lewis SR. QTL mapping of morphological characteristics that correlated to drought tolerance in St. Augustinegrass. PLOS ONE, 2022;17:e0268004.
Schoonmaker AN. Analysis and computational methods for understanding disease resistance at the genomic level in complex plant systems [dissertation]. Raleigh (US): North Carolina State University; 2023.
Chen L, Fan J, Hu L, Hu Z, Xie Y, Zhang Y, et al. A transcriptomic analysis of bermudagrass (Cynodon dactylon) provides novel insights into the basis of low temperature tolerance. BMC Plant Biol., 2015;15:216.
Shu YJ, Zhou YL, Mu KB, Hu HM, Chen M, He QY, et al. A transcriptomic analysis reveals soybean seed pre-harvest deterioration resistance pathways under high temperature and humidity stress. Genome, 2020;63:115-24.
Brown JM, Weldt CE, Holloway HMP, Tuong TD, Patton AJ, DaCosta M, et al. Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation. Grass Research, 2023;3:25.
Gelli M, Konda AR, Liu K, Zhang C, Clemente TE, Holding DR, et al. Validation of QTL mapping and transcriptome profiling for identification of candidate genes associated with nitrogen stress tolerance in sorghum. BMC Plant Biol, 2017;17:123.
Diouf I, Albert E, Duboscq R, Santoni S, Bitton F, Gricourt J, Causse M. Integration of QTL, transcriptome and polymorphism studies reveals candidate genes for water stress response in tomato. Genes, 2020;11:900.
Blum A. Plant Breeding for Water-Limited Environments. London: Springer; 2010. p. 1-210.
Holland JB. Genetic architecture of complex traits in plants. Curr. Opin. Plant Biol. 2007;10:156-61.
Yu X, Mulkey SE, Zuleta MC, Arellano C, Ma B, Milla-Lewis SR. Quantitative trait loci associated with gray leaf spot resistance in St. Augustinegrass. Plant Dis. 2020;104:2799-806.
Madden LV, Hughes G, Van Den Bosch F. The study of plant disease epidemics. St. Paul: APS Press; 2007.
Zhang J, Virk S, Porter W, Kenworthy K, Sullivan D, Schwartz B. Applications of unmanned aerial vehicle based imagery in turfgrass field trials. Frontiers in Plant Science, 2019;10:279.
Rockstad GB, Austin RE, Gouveia BT, Carbajal EM, Milla‐Lewis SR. Assessing unmanned aerial vehicle‐based imagery for breeding applications in St. Augustinegrass under drought and non‐drought conditions. Crop Science, 2024;64:496-510.
Merewitz E, Belanger F, Warnke S, Huang B, Bonos S. Quantitative trait loci associated with drought tolerance in creeping bentgrass. Crop Science, 2014;54:2314-24.
Jespersen D, Merewitz E, Xu Y, Honig J, Bonos S, Meyer W, Huang B. Quantitative trait loci associated with physiological traits for heat tolerance in creeping bentgrass. Crop Science, 201;56:1314-29.
Zhang D, Zhang H, Chu S, Li H, Chi Y, Triebwasser-Freese D, et al. Integrating QTL mapping and transcriptomics identifies candidate genes underlying QTLs associated with soybean tolerance to low-phosphorus stress. Plant Molecular Biology, 2017;93:137-50.
Ye J, Yang Y, Chen B, Shi J, Luo M, Zhan J, et al. An integrated analysis of QTL mapping and RNA sequencing provides further insights and promising candidates for pod number variation in rapeseed (Brassica napus L.). BMC Genomics, 2017;18:71.
Zhang Q, Liu H, Wu X, Wang W. Identification of drought tolerant mechanisms in a drought-tolerant maize mutant based on physiological, biochemical and transcriptomic analyses. BMC Plant Biology, 2020;20:315.
Yennamalli RM, Rader AJ, Kenny AJ, Wolt JD, Sen TZ. Endoglucanases: Insights into thermostability for biofuel applications. Biotechnology for Biofuels, 2013;6:136.
Muthurajan R, Rahman H, Manoharan M, Ramanathan V, Nallathambi J. Drought responsive transcriptome profiling in roots of contrasting rice genotypes. Indian Journal of Plant Physiology, 2018;23:393-407.
Xu C, Huang B. Comparative analysis of drought responsive proteins in Kentucky bluegrass cultivars contrasting in drought tolerance. Crop Science, 2010;50:2543-52.
Halušková L, Valentovičová K, Huttová J, Mistrík I, Tamás L. Effect of abiotic stresses on glutathione peroxidase and glutathione S-transferase activity in barley root tips. Plant Physiology and Biochemistry, 2009;47:1069-74.
Park WJ, Campbell BT. Aquaporins as targets for stress tolerance in plants: Genomic complexity and perspectives. Turkish Journal of Botany, 2015;39:879-86.
Ding YF, Cao JM, Ni L, Zhu Y, Zhang AY, Tan MP, Jiang MY. ZmCPK11 is involved in abscisic acid-induced antioxidant defence and functions upstream of ZmMPK5 in abscisic acid signalling in maize. J. Exp. Bot. 2013;64:871-84.
Kumar K, Rao KP, Sharma P, Sinha AK. Differential regulation of rice mitogen activated protein kinase kinase (MKK) by abiotic stress. Plant Physiol. Biochem. 2008;46:891-97.
Wang M, Yue H, Feng KW, Deng PC, Song WN, Nie XJ. Genome-wide identification, phylogeny and expressional profiles of mitogen activated protein kinase kinase kinase (MAPKKK) gene family in bread wheat (Triticum aestivum L.) BMC Genom. 2016;17:668.
Zhang XY, Wang LM, Xu XY, Cai CP, Guo WZ. Genome-wide identification of mitogen-activated protein kinase gene family in Gossypium raimondii and the function of their corresponding orthologs in tetraploid cultivated cotton. BMC Plant Biol. 2014;14:345.
Liu Y, Zhou MY, Gao ZX, Ren W, Yang FL, He H, Zhao JR. RNA-seq analysis reveals MAPKKK family members related to drought tolerance in maize. PLoS ONE. 2015;10:e0143128.
Krasensky J, Jonak C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. Journal of Experimental Botany, 2012;63:1593-608.
Thalmann M, Santelia D. Starch as a determinant of plant fitness under abiotic stress. New Phytologist, 2017;214:943-51.
Lu P, Magwanga RO, Kirungu JN, Hu Y, Dong Q, Cai X, et al. Overexpression of cotton a DTX/MATE gene enhances drought, salt, and cold stress tolerance in transgenic Arabidopsis. Front Plant Sci. 2019;10:299.
Tran LSP, Shinozaki K, Yamaguchi‐Shinozaki K. Role of cytokinin responsive two‐component system in ABA and osmotic stress signaling. Plant Signal. Behav. 2010;5:148-50.
Wang R, Liu L, Kong Z, Li S, Lu L, Kabir N, Chen G, et al. Identification of GhLOG gene family revealed that GhLOG3 is involved in regulating salinity tolerance in cotton (Gossypium hirsutum L.). Plant Physiol. Biochem. 2021;166:328-40.
de Mendiburu F. agricolae: Statistical procedures for agricultural research. R package version 1.3-5. 2021. https://CRAN.R-project.org/package=agricolae.
Wei T, Simko V. R package 'corrplot': Visualization of a correlation matrix. R package version 0.92. 2021. https://github.com/taiyun/corrplot.
Butler DG, Cullis BR, Gilmour AR, Gogel BG, Thompson R. ASReml-R reference manual version 4.2. 2018. https://asreml.kb.vsni.co.uk/wp-content/uploads/sites/3/ASReml-R-Reference-Manual-4.2.pdf.
Afanador L, Haley SD, Kelly JD. Adoption of a mini-prep DNA extraction method for RAPD marker analysis in common bean. Bean Improv Coop, 1993;35:10-11.
Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, Mitchell SE. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLOS ONE, 2011;6:e19379.
Melo ATO, Bartaula R, Hale I. GBS-SNP-CROP: A reference-optional pipeline for SNP discovery and plant germplasm characterization using variable length, paired-end genotyping-by-sequencing data. BMC Bioinformatics, 2016;17:29.
Chen S, Zhou Y, Chen Y, Gu J. fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 2018;34:i884-90.
Mollinari M, Garcia AAF. Linkage analysis and haplotype phasing in experimental autopolyploid populations with high ploidy level using hidden Markov models. G3 Genes|Genomes|Genetics, 2019;9:3297-314.
Preedy KF, Hackett CA. A rapid marker ordering approach for high-density genetic linkage maps in experimental autotetraploid populations using multidimensional scaling. Theoretical and Applied Genetics, 2016;129:2117-32.
da Silva Pereira G, Gemenet DC, Mollinari M, Olukolu BA, Wood JC, Diaz F, et al. Multiple QTL mapping in autopolyploids: a random-effect model approach with application in a hexaploid sweetpotato full-sib population. Genetics, 2020;215:579-95.
Weldt MC. From transcriptomics to QTL mapping: exploring drought resistance in St. Augustinegrass [thesis]. Raleigh (US): North Carolina State University; 2024.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494-512.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139-40.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology, 2010;11:R14.
Voorrips RE. MapChart: Software for the graphical presentation of linkage maps and QTLs. The Journal of Heredity, 2002;93:77-8.

No competing interests reported.

Additionalfile1.pdf
Additional file 1.pdf: Significance tests for each trait and relevant variance components across years.
Additionalfile2.pdf
Additional file 2.pdf: Synteny between genomic marker order and linkage group marker order.
Additionalfile3.xlsx
Additional file 3.xlsx: List of differentially expressed genes (DEGs) in St. Augustinegrass in response to water stress.
Additionalfile4.xlsx
Additional file 4.xlsx: Gene Ontology (GO) classification of differentially expressed genes (DEGs) in St. Augustinegrass in response to water stress.
Additionalfile5.xlsx
Additional file 5.xlsx: KEGG pathway enrichment of differentially expressed genes (DEGs) in St. Augustinegrass in response to water stress.

Download PDF

Reviews received at journal
31 Aug, 2024
Reviewers agreed at journal
24 Aug, 2024
Reviewers agreed at journal
24 Aug, 2024
Reviewers agreed at journal
23 Aug, 2024
Reviewers invited by journal
15 Jul, 2024
Editor assigned by journal
15 Jul, 2024
Submission checks completed at journal
15 Jul, 2024
First submitted to journal
07 Jul, 2024

You are reading this latest preprint version

Integration of QTL and Transcriptome Studies Reveals Candidate Genes for Water Stress Response in St. Augustinegrass

Status:

Version 1

Abstract

Background

Results

Conclusions

Figures

BACKGROUND

RESULTS

Drought-related traits evaluation

SNP discovery and linkage map construction

QTL detection for drought-related traits

QTL overlap across studies

Transcriptome assembly

Differentially expressed genes (DEGs) under water stress

GO classification of differentially expressed genes

KEGG pathway enrichment of differentially expressed genes

Co-localization of the DEGs and drought QTL

DISCUSSION

CONCLUSION

MATERIALS AND METHODS

Mapping population development

Greenhouse drought evaluations

Statistical analysis

Library construction and sequencing

SNP calling

Linkage map construction

QTL mapping

Materials and treatment for transcriptomics study

RNA extraction, library preparation and sequencing

Quality control, transcriptome assembly and gene annotation

Differential expression analysis

GO and KEGG enrichment analysis

QTL validation with previous studies

Co-Localization of the DEGs and drought QTL

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1