Quantitative trait loci mapping and candidate gene analysis of stoma-related traits in wheat (Triticum aestivum L.) glumes

Test material

A wheat double haploid (DH) population (Liu et al., 2013), including 150 lines that derived from a cross between Hanxuan 10 and Lumai 14 was used in this study. All the 150 lines and parents were grown at the experimental farm (37°25′N, 112°35′E) of Shanxi Agricultural University in 2018 and 2019. The field experiments were conducted by randomized complete block design (RCBD) with three replicates. Each plot consisted of two rows of 2 m in length, with 0.25 m between rows. Water and fertilizer management during the growth period was complied with the local production practice.

Measurement and calculation of stomata-related traits in glumes

For each DH line and parents, three flowering plants with consistent growth were tagged. The middle spikes of the three plants were quickly placed into a 2-mL centrifuge tube with FAA solution (formalin: acetic acid: 70% alcohol = 1:1:18) three days after anthesis, respectively. The sampling was started at 10 am on the third day after anthesis. The average temperature at the time of sampling was 13 and 17 °C, and the moisture content of topsoil (0∼30 cm) was approximately 15% and 13% (v/v) in 2018 and 2019, respectively. Stomatal density (SD), stomatal length (SL), and stomatal width (SW) on the top, middle, and base of wheat glumes were measured as previously described in Wang et al. (2018).

Calculations of other stoma-related traits were as follows:

Potential conductance index (PCI) = SL² × SD × 10⁻⁴ (Nicholas & Richardson, 2009)

Stomatal area (SA) = π × 1/2 SL × 1/2 SW (Robert, Gretchen & Richard, 2000)

Stomatal relative area (SRA) = SD × SA × 10⁻⁴ (Sun et al., 2021)

Data analysis

The relevant t-test and analysis of variance (ANOVA) were carried out by the statistical software package SPSS v.17.0; the frequency distribution map was generated by Excel 2007; and the phenotypic correlation analysis map was created in the R-package corrplot (Wei & Simko, 2013).

QTL mapping

The genetic map of the DH population was constructed by Jing Ruilian’s team at the Institute of Crop Science, Chinese Academy of Agricultural Sciences (Li et al., 2019). The linkage map is 4,082.4 cM in length and contains 1630 SNP and 224 SSR markers, with 2.2 cM per bin on average (Shi et al., 2020). QTL mapping was performed as previously described in Li et al. (2021) using the IciMapping 4.1 software. The LOD score for declaring a QTL was 2.5 for each trait.

Prediction of candidate genes

Candidate genes in associated loci were predicted according to the reference genome sequence of ‘Chinese Spring’ wheat (IWGSC RefSeqv1.1) published by the International Wheat Genome Sequencing Consortium. Gene annotation was carried out by referring to the Ensembl Plants database (https://plants.ensembl.org/index.html). We used a publicly available database, WheatOmics (http://wheatomics.sdau.edu.cn/) (Ma et al., 2021), to obtain the expression profiles of all candidate genes.

Phenotypic variation of glume stoma-related traits in the DH population and parents

We observed regular rows of stomata on the top, middle, and base of glumes in parents (Fig. 1). SD of different parts of glumes in Lumai 14 showed significant or highly significant differences between the two years, while Hanxuan 10 had a significant difference in SD only at the base of glumes (Table 1). For DH lines, except for SD at the base of glumes, the phenotypic values of SD of other parts of glumes in 2019 were significantly lower than those in 2018 (Table 1). These results indicated that SD of glumes is greatly affected by the environment, and this phenomenon is more obvious in Lumai 14 than Hanxuan 10. In general, from the top to the base, SD of glumes gradually became less (Table 1).

Except for the middle of glumes in Hanxuan10, differences in SL of the rest parts of glumes between the two years were not significant in the two parents. Except for the top of glumes of Lumai 14, differences in SW of the rest parts of glumes between the two years were not significant in the two parents. For the DH lines, the difference in SL of glumes was not significant between the two years, but SW in 2019 was highly significantly smaller than that in 2018 (Table 1). These results indicated that SW of glumes is more affected by the environment than SL. For the 150 DH lines, there were genetic differences in the same stoma-related traits among different lines, and this difference varied among different traits. The coefficient of variation of SD in the DH population was greater than 10 in both 2018 and 2019, while the coefficient of variation of both SL and SW was less than 10. This also showed that different stomatal-related traits are affected differently by the environment. In addition, SL and SW did not change significantly from the top to the base of wheat glumes.

Except for the middle of glumes of Lumai 14, differences in stomatal PCI of the rest parts of glumes between the two years were not significant in the two parents. Differences in stomatal PCI of all parts between the two years were not significant for the DH lines neither. Differences in SA of all parts between the two years were not significant in the two parents. Except for the top and middle of glumes of Lumai 14, differences in SRA of the rest parts of glumes between the two years were not significant in the two parents (Table 1).

Stoma-related traits of the DH lines showed continuous transgressive segregation with skewness and kurtosis values close to zero, suggesting normal distribution. All target traits were thus quantitatively controlled by multiple genes and were suitable for QTL mapping (Table 1) (Fig. S1).

Correlation between glume stoma-related traits

For all parts of glumes, SD showed a highly significantly negative correlation with SL in 2018 and 2019; however, the correlation between SD and SW was not significant. There was a significantly positive correlation between SL and SW in 2018, but such correlation was not significant in 2019. In addition, SD showed a highly significant and positive correlation with PCI and SRA in 2018 and 2019. SL showed a significantly positive correlation with SA in both years. SW showed a significantly positive correlation with SA and SRA in both years. In 2018, any two of PCI, SA, and SRA were significantly positively correlated, while the degree of correlation was weakened in 2019 (Fig. 2 Table S1).

QTL mapping for traits related to stomata in wheat glumes

A total of 61 QTLs for traits related to stomata of wheat glumes were detected in the two years. The phenotypic variation of these QTLs ranged from 3.63 to 19.02%. The LOD score ranged from 2.51 to 22.12, and the QTLs were distributed on 16 chromosomes including 1A, 1D, 2A, 2B 2D, 3A, 3D, 4B, 5A, 5B, 5D, 6A, 6B, 7A, 7B, and 7D, respectively (Fig. 3 Table 2). A total of nine QTLs were detected in both years.

A total of 10 QTLs corresponding to SD of glumes were detected in the two years. The phenotypic variation of these QTLs ranged from 6.78 to 11.41%. Among these QTLs, QSDt-2A was detected in both years; QSDt-5A and QSDm-5A were detected in the same interval in 2018. Twenty-one QTLs were associated with stomatal size and detected in both years. Among these QTLs, QSLt-2D, QSLm-7A, and QSLb-6A-1 were detected in both years. QSLb-6A-2 and QSWb-6A were detected in the same interval in 2019. A total of 30 QTLs associated with PCI, SA, and SRA were detected in the two years. Among these QTLs, QPCIb-6A-1, QSAt-2D, QSAm-3A, and QSRAb-6A-2 were detected in the two years. Among all QTLs, QSRAb-6A-1 had the largest LOD value (22.12) (Fig. 3 Table 2).

Among all QTLs, two QTL hotspots were found on chromosome 6A. One was in the interval AX-109400932∼AX-110985652, which contained four QTLs related to SL, SD, PCI, and SRA in the two years. The other was in the interval AX-108972184∼AX-108752564, which contained four QTLs related to SW, SD, PCI, and SRA in the two years (Fig. 3 Table 2).

The prediction of candidate genes

For the first QTL hotspot, the physical positions of the two markers AX-109400932 and AX-110985652 were 73571398 and 76990896 bp, respectively. According to the reference genome sequence of ‘Chinese Spring’ wheat (IWGSC RefSeqv1.1), a total of 33 genes were found between the two markers, and gene annotation was carried out by referring to the Ensembl Plants database (Table S2). For the other QTL hotspot, the physical positions of the two markers AX-108972184 and AX-108752564 were 83931623 and 86272494 bp, respectively, covering a total of 24 genes (Table S2). Using WheatOmics, the expression levels of these 57 genes in wheat glumes and spikes were predicted (The International Wheat Genome Sequencing Consortium IWGSC, 2018). The results showed that there were five genes (TraesCS6A01G105400, TraesCS6A01G106400, TraesCS6A01G115100, TraesCS6A01G115400, and TraesCS6A01G116200) with higher expression levels in wheat glumes or spikes (Fig. 4 Table 3).

Then, the expression levels of these five genes under different abiotic stresses were predicted (Oono et al., 2013; Liu et al., 2015; Nazanin et al., 2019). The results showed that the expression of these genes was induced by different abiotic stresses, and the expression patterns of different genes under the same abiotic stress were also different (Fig. 5). For example, the expression of TraesCS6A01G105400 was significantly reduced after being subjected to low phosphorus stress in wheat shoots and roots. The expression levels of TraesCS6A01G115100 and TraesCS6A01G116200 both increased after being subjected to low phosphorus stress in shoots and roots. The expression of TraesCS6A01G106400 was only significantly increased in the root after being subjected to low phosphorus stress. The expression level of TraesCS6A01G105400 was significantly reduced in wheat seedling leaves subjected to drought stress, heat stress, and their combination. The expression levels of TraesCS6A01G106400 and TraesCS6A01G116200 were significantly increased after six hours of drought stress, heat stress, and their combination. Among these five genes, only the expression level of TraesCS6A01G115400 was significantly decreased after being exposed to salt stress (Fig. 5).

Phenotypic correlation of stoma-related traits in wheat glumes

Stomata are the channel for water and gas exchange in the process of wheat photosynthesis and respiration, which indirectly affect the yield of wheat (Berger & Altmann, 2000). Studies have shown that a variety of environmental factors at different growth and developmental stages of plants can affect the formation of stomata, such as water (Stephens & Waugh, 2017), temperature (Qi & Torii, 2018), light (Boccalandro et al., 2009), and CO₂ concentration (Hu et al., 2010). Aasamaa, Sober & Rahi (2001) reported that the stomatal length of forest tree species decreased with increasing drought. For some light-loving crops, the formation of stomata can be promoted by increasing the light intensity. In the present study, stomatal density and stomatal width of wheat glumes showed significant differences between the two years, but stomatal length showed no significant differences (Table 1). These results suggest that stomatal length may have higher stability in response to different environmental conditions than stomatal density and width. In addition, we also compared the differences of the same traits in different parts of wheat glumes, and the results showed that from the top to the base of glumes, the stomatal density gradually decreased, but the stomatal length and width did not change significantly. Moreover, four of the nine stable QTLs mapped were associated with stoma-related traits at the base of wheat glumes. It is speculated that this may be because the stomatal properties at the base of the glume are the most stable compared with other parts.

Franks & Beerling (2009) have reported that the negative correlation between stomatal size and stomatal density helps to adjust the plasticity of stomata, thereby regulating the maximal stomatal conductance of wheat. There was a significantly negative correlation between stomatal density and stomatal size in wheat leaves (Wang et al., 2016). In addition, similar phenomena have been observed in other crops. (Ishimaru et al., 2001; Ohsumi et al., 2007). The present results showed that stomatal density was negatively correlated with stomatal length in each part of wheat glumes in the two years. However, there was no significant correlation between stomatal density and stomatal width. Previous studies found that stomatal length and stomatal width in wheat leaves were significantly positively correlated (Wang et al., 2016). The present results showed that stomatal length and stomatal width were significantly positively correlated in 2018 (Fig. 2 Table S1). Therefore, wheat glumes can improve their adaptability to different environmental conditions by coordinating the relationship among stomatal density, stomatal length, and stomatal width.

Pleiotropy of QTLs for stoma-related traits in wheat

Various studies have found that QTLs of closely related traits may be located on the same or nearby positions on the chromosomes (Fracheboud et al., 2002; Tuberosa et al., 2002). In the present study, QSLb-6A-1 corresponding to stomatal length at the base of glumes, QSDb-6A for stomatal density at the base of glumes, QPCIb-6A-1 for stomatal PCI at the base of glumes, and QSRAd-6A-1 for stomatal relative area at the base of glumes were detected in the interval of AX-109400932∼AX-110985652 on chromosome 6A in 2018 and 2019. In the vicinity of this interval, QSWb-6A, QSLb-6A-2, QPCIb-6A-2, and QSRAb-6A-2 were detected in the interval AX-108972184∼AX-108752564 on chromosome 6A in the two years. QSDt-5A, QSDm-5A, and QSRAt-5A were all located in the interval AX-111662464∼AX-95683796. QSDt-5D and QSLt-5D were located in the interval Xgwm205.2∼AX-89390905. QSDm-1A and QSLt-1A were located in the interval AX-111105973∼AX-94402739 (Fig. 3 Table 2).

In addition, compared with previous studies, we found a QTL QSWt-4B in this study was located within the region of QAGsw4B in the previous study (Wang et al., 2018). QSAb-5A-2 was located within the region of QAGsd5A, QMGsd5A-2, and QAGsl5A. (Wang et al., 2018). QSRAm-5B was located close to the region of QPsd5B, and QSDd-5A was located within the region of QSD5A-2 (Wang et al., 2016).

Therefore, the QTLs for above-mentioned stoma-related traits, which were detected in different parts of wheat, different growth periods, and various environments, were significant markers for stoma-related traits in wheat. Furthermore, these findings implied stomatal density and size of wheat leaves and glumes may be controlled by the same or pleiotropic genes. The markers that were localized within a QTL interval associated with stoma-related traits not only validated the QTL but also provided more closely linked markers. These markers will be useful to reveal advanced wheat varieties in wheat breeding programs based on marker-assisted selection approaches.

Prediction of candidate genes related to stomata in wheat glumes

In this study, 57 genes were found in two intervals, AX-109400932∼AX-110985652 (physical range 73571398-76990896 bp) and AX-108972184∼AX-108752564 (physical range 83931623-86272494 bp), on chromosome 6A, and five candidate genes (i.e., TraesCS6A01G105400, TraesCS6A01G106400, TraesCS6A02G115100, TraesCS6A02G115400, and TraesCS6A02G116200) were screened out, according to their expression levels in wheat glumes or spikes (Fig. 4 Table 3).

The expression level of TraesCS6A01G105400 in wheat glumes and spikes was the highest among all candidate genes (Fig. 4), and its functional annotation was 50S ribosomal protein L3 (Table 3). The TraesCS6A01G105400 expression decreased significantly under low phosphorus, drought, and heat stress conditions (Fig. 5). The functional annotation of the candidate gene TraesCS6A01G106400 is stress-associated endoplasmic reticulum (ER) protein (Table 3). The ER plays a crucial role in the maintenance of cellular homeostasis. ER stress is a widely existed stress mechanism to external stimuli in plants and animals. This pathway maintains the ER homeostasis and alleviates stress damage through regulation of a series of gene expressions (Park & Park, 2019). The expression of TraesCS6A01G106400 was significantly increased under salt stress (Fig. 5). The functional annotation of the candidate gene TraesCS6A01G115100 is purple acid phosphatase (PAP) (Table 3), and its expression was significantly increased under low phosphorus stress (Fig. 5). PAPs are members of the metallo-phosphoesterase family identified from a wide range of plants. PAPs have mostly been studied for their potential involvement in phosphorus acquisition and redistribution because of their ability to catalyze the hydrolysis of activated phosphate esters and anhydrides under acidic conditions (Olczak, Morawiecka & Watorek, 2003). Recent studies also showed that PAPs play important roles in modulating plant carbon metabolism, cell wall synthesis, pathogen resistance, etc (Kaida et al., 2009; Sun et al., 2012; Zhang et al., 2014). TraesCS6A01G115400 was specifically expressed in wheat glumes, and its functional annotation is calcium-dependent lipid-binding family protein (Table 3). Ca²⁺ is a secondary messenger in plants that regulates virtually all aspects of plant development and responses to environmental stimuli. Ca²⁺ tends to rapidly rise under abiotic stresses (Bartels & Sunkar, 2005). Several proteins have been reported to be activated or translocated in the presence of Ca²⁺ including Calcium-dependent lipid-binding protein (Hurley & Misra, 2000). The expression of TRAESCS6A01G115400 in roots of wheat seedlings was significantly decreased under salt stress (Fig. 5). The functional annotation of the candidate gene TraesCS6A01G116200 is ATP-dependent RNA helicase (Table 3). ATP-dependent RNA helicase can be found in many organisms, which is involved in the multi-dimensional metabolism of RNA and plays an important role in plant growth and development, especially in abiotic stress response (Kim et al., 2008). The expression of TraesCS6A01G116200 was significantly increased under low phosphorus, drought, and heat stress (Fig. 5).

The above five candidate genes not only had high expression levels in wheat glumes or spikes, but are also induced by a variety of abiotic stresses. Studies have shown that when plants are subjected to abiotic stress, they can change their photosynthetic rate and transpiration rate by adjusting the stomata size, stomata density, and stomata distribution to deal with bad external environments (Doheny-Adams et al., 2012; Franks et al., 2015). Therefore, it will be greatly helpful to analyze whether these five genes affect the formation of wheat stomata through transgenic experiments in future research, and to explore the mechanism of their functions.