Identification of genetic loci for early maturity in spring bread wheat using the association analysis and gene dissection

Plant material and phenotyping

In this study a collection of 92 Russian spring wheat varieties was investigated that originated from the breeding centers of seven regions such as Samara, Altai, Novosibirsk, Krasnoyarsk, Kemerovo, Tyumen, and Omsk. The varieties were maintained in the GeneAgro seed bank of Institute of Cytology and Genetics of Siberian Branch of the Russian Academy of Sciences (ICG SB RAS) (https://ckp.icgen.ru/plants/). The wheat varieties had been adapted to the climatic conditions of Western Siberia whose summer period is rather short, so different environmental troubles often happen at the end of the vegetation period. More detailed information about the wheat varieties’ origin is presented in Leonova et al. (2022).

For mapping the genetic loci, we obtained two phenotyping datasets for this panel of spring wheat varieties. The first dataset was obtained in 2018 and described the heading, heading to maturity and maturity periods to cover the total seedling to maturity time. The second dataset included 10-year data (from 2005 to 2014) evaluating the time from seedling to maturity.

The plants were cultivated in the two experimental fields of ICG SB RAS in the Novosibirsk Region (field 1: 54.884987°N, 82.899519°E; field 2: 54.914070°N, 82.975379°E). The accessions were phenotyped on both fields in 2018 and on field 1 in 2005–2014.

The phenological stages were evaluated following the Zadoks classification (Zadoks, Chang & Konzak, 1974). The heading time was recorded at the moment when 70% of the plants in the plot had emerged with the ear from the tube by one-half, and it was calculated as the number of days from seedling to heading (corresponding to stage GS55). The maturity was determined by hard grains, yellowing, and drying of the ears and stems (corresponding to stage GS92).

The seeds were sown in the second half of May in 2 repetitions on plots of 2 m². The soil in both fields was leached chernozem. The humus thickness varied from 40 to 60 cm, and the humus content was 4.2%. The soil was slightly acidic (pH 6.7–6.8), with its nitrogen content being 0.34%, total phosphorus-0.30% and potassium-0.13%.

To test the effect of identified SNPs, an F_2:3 population was obtained from a cross of mid-season cultivar Obskaya 2 and early-maturing cultivar Novosibirskaya 15. The population was phenotyped in 2018 in field 1 in 2 repetitions on plots of 0.25 m².

The weather data (average daily temperature and humidity by months) provided by a meteorological station located in the Novosibirsk region at 54.90 N, 82.95 E and an altitude of 131 m were taken from http://www.pogodaiklimat.ru. The air temperature and precipitation amount during the experiments were contrasting both during one growing season and over the years. The most optimal weather to ensure high productivity was registered in 2005 and 2009, whose temperatures and precipitation amounts in these years were closest to the mean seasonal values. The moisture supply was insufficient in 2008 (166.8 mm), 2010 (135.9 mm), 2011 (153.5 mm), 2012 (102.7 mm) and 2014 (176 mm). An extremely severe drought was noted in July 2012, with 3.7 mm of rain falling during the month, which was only 6.1% of the mean seasonal value.

Growing degree days (GDD) from seedling to maturity, were calculated according to the formula of McMaster (1997b): GDD = (Tmax + Tmin)/2 – Tbase, where Tmax and Tmin are the maximum and minimum temperatures for the day and Tbase was 0 °C (McMaster, 1997a). The GDD was highest in 2012 (more than 2,200). Cold and rainy weather at the beginning of the growing season was observed in 2007, 2010, 2013, 2014 and 2018 with GDD less than 1,600. The most optimal air temperature for plant development was in 2008 and 2014 when the GDD varied from 1,800 °C to 1,900 °C. The weather was especially unfavourable for seedlings in 2018 when the average air temperature was 4 °С below the norm in May and the precipitation exceeded the mean seasonal value by 126.8%.

DNA extraction and high-throughput genotyping

The DNA was extracted from young leaves using the modified sodium bisulfite protocol described by Kiseleva et al. (2016) and purified using the Bio-Silica Kit for DNA Purification from Reaction Mixtures as per the manufacturer’s protocol. The samples were genotyped using the Illumina Infinium 15 K Wheat containing 13,007 SNPs by TraitGenetics GmbH (www.traitgenetics.com).

Statistical analysis

Analysis of variance (ANOVA) to explore the association between traits was performed using the ‘aov’ function in R. Spearman’s correlation coefficients (r²) were calculated. The model considering genotypic effects as fixed and environmental effects as random was used to calculate the best linear unbiased estimates (BLUEs) of each genotype. To check if the distribution of the traits was normal, we used the Shapiro–Wilk normality test (R function). To estimate the maturity and heading time heritability, all effects except the intercept were assumed as random, so the component variance was estimated using R package lme4 (Bates et al., 2020) and H² was calculated using the “Standard” broad-sense heritability method (Holland, Nyquist & Cervantes-Martínez, 2003).

Association analysis

Only polymorphic markers (SNPs) were selected for subsequent analysis. The genotyping data were filtered by minor allele frequency (MAF > 0.05) and by missing data less than 5% of population. Population structure analysis was performed using the LEA algorithm (snmf function with K = 1–10) (Frichot & François, 2015).

A mixed linear model (MLM) was performed in TASSEL V5.1 to carry out the association analysis, which accounted for fixed effects of SNPs and population structure (Q) and random effects for kinship (K). Quantile-quantile and Manhattan plots were generated using the R package “GWASTools” (Conomos et al., 2012).

Linkage disequilibrium (LD) between SNP markers was calculated using the R package “genetics” (Warnes et al., 2019). LD decay plots were generated using R package “LDheatmap” (Shin et al., 2006). The chromosome positions of SNP markers were established according to the IWGSC Reference Sequence v1.0 Annotation (Consortium IWGS, 2018).

Annotation of candidate genes

To provide insight into the functional specification of candidate genes, they were annotated using the IWGSC 1.0 data. To assess the gene expression, the developmental time course of the common wheat cultivar Azhurnaya (Ramírez-González et al., 2018) deposited in Wheat Expression Browser expVIP (Borrill, Ramirez-Gonzalez & Uauy, 2016) was used. Visualization and clustering were performed using R package d3heatmap. Gene-specific primers were designed in PrimerQuest Tool software (Integrated DNA Technologies, Inc., Coralville, IA, USA) (Table S1).

RNA isolation and RT-qPCR analysis

According to the results of the association study, 10 wheat varieties with different haplotypes in the QTL located on 7B chromosome were selected. The varieties considered as early and late maturing and having the respective CCGACTGACT and TAAGTGAGCT haplotypes were cultivated in the year 2020 on the experimental field of ICG SB RAS in Novosibirsk. The plants were harvested at the start of the anthesis stage GS61 (Zadoks, Chang & Konzak, 1974). The stage was chosen based of the candidate gene expression data from Wheat Expression Browser expVIP (Borrill, Ramirez-Gonzalez & Uauy, 2016). Three replicates from each genotype were harvested into liquid nitrogen in 7 h after the beginning of the light period (about 10 a.m.). Each biological replicate consisted of three plants of each cultivar combined into a bulk, with each combination replicated three times (biological replicates) for a total of 30 plants sampled.

The RNA was extracted using Plant RNA MiniPrep (Zymo Research, Tustin, CA, USA), followed by DNase treatment with a RNase-Free DNase set (QIAGEN, Hilden, Germany). The quantity and quality of the extracted RNA were assessed using a QUBIT 4 fluorometer with RNA BR and RNA IQ kits respectively (Thermo Fisher Scientific, Vilnius, Lithuania). cDNA was synthesized using RevertAid First Strand cDNA Synthesis (Thermo Fisher Scientific, Vilnius, Lithuania) following the manufacturer’s protocol with 2 μg of total RNA as a template and random hexamer primers (Thermo Fisher Scientific, Vilnius, Lithuania). In total, 2 μl of 20-fold final cDNA dilution was used for the following analysis.

The fluorescence data were collected using a QuantStudio™ 5 Real-Time PCR System (Thermo Fisher Scientific, Vilnius, Lithuania) with PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, Vilnius, Lithuania). The measurements were performed in three technical replicates. The reaction products were checked using melting curve analysis.

To detect the most applicable reference genes, the expression of seven candidate reference genes selected according to the literature (Table S2) was analyzed. Applying the geNorm algorithm (Vandesompele et al., 2002) integared in qBase+ software (Hellemans et al., 2008) two genes (MetAP1 and SAR) were chosen as the best references (geNorm V < 0.15).

The PCR efficiencies were determined using the LinReg software (Ruijter, Ilgun & Gunst, 2014). The Calibrated Normalized Relative quantities (CNRQ) were calculated in qBase+ (Hellemans et al., 2008). The target-gene expression was normalized against two reference genes (MetAP1 and SAR). The plots representing expression ratios were generated with the ggplot2 (Wickham, 2016) and ggpubr (Kassambara, 2023) R packages.

Phenotyping

The 10-year period data included information only on the days from seedling to maturity (maturity time). To reveal correlations between the heading and maturity phenological phases, the times from seedling to heading (heading time, HT), from heading to maturity (HTMT), and the total time from seedling to maturity (maturity time, MT) were evaluated for two environmental conditions in 2018. Using these data, HT and HTMT contributions to MT were estimated.

The results of this experiment demonstrated that the time from seedling to heading (HT) of the studied varieties ranged from 35 to 44 days, with a heritability of 53% for this trait. While the total time from seedling to maturity (MT) varied from 77 to 104 days (Fig. 1). The large MT variation and extension (about a 4-week difference) may be due to unfavorable weather conditions in 2018, when the temperature in May was 4 °С below the norm, and the precipitation during growing season significantly exceeded the average.

Correlation analysis of these parameters demonstrated that HTMT strongly correlated with MT (r = 0.96, p-value < 2.2E−16), while HT moderately correlated with MT (r = 0.58, p-value = 2.089E−10) (Fig. 1).

Ten-year evaluation (2005–2014) of MT duration showed that it varied from 60 to 96 days in different years (Table S3). ANOVA indicated that such parameters as the genetic component and environment both significantly contributed to MT (Table S4). Individual variance component analysis of the ten-year data demonstrated that the genotypic variance principally contributed to MT (73.5%).

Five varieties (Polushko, Novosibirskaya-15, Salimovka, Novosibirskaya 22, Lutescense 25) had the earliest maturity time (MT < 72 days) if compared with cv. Novosibirskaya 81, Omskaya 28, Tulaikovskaya zolotistaya, Tulaikovskaya 10, and Baganskaya 93 which matured later than the others (MT > 81 days). These varieties were subsequently used to study the expression of candidate genes as contrasting genotypes.

For further analysis and minimization in assessing the environmental impact, the best linear unbiased estimates based on the data for 11 years of observation were estimated (2005–2014, 2018).

Association analysis

Association mapping for MT was performed using the BLUEs of 11-year phenotype estimation of 96 wheat varieties and its results are summarized in Table 1. A total of 22 significant SNPs were located on chromosomes 2A, 3B, 4A, 5B, 7A, and 7B (Fig. 2A, Table 1). Most of these chromosomes are represented by one or two significant markers except for chromosome 7BL where ten SNPs were tightly linked to each other (Fig. 2C, Table 1). In the following analysis, only the 7BL locus was considered.

LD analysis based on the r² measure between significant SNPs (p < 0.001) demonstrated that tightly linked markers within a QTL on the long arm of chromosome 7B were located within a 658,219,587–680,161,741 bp interval on the physical map. Excalibur_c31707_302 was distant from other SNPs (about 21 Mb), however, it was found to be tightly linked genetically (r² ~ 0.9–1). According to GWAS, the 7B locus negatively influenced the trait’s manifestation by reducing MT, whereas markers Kukri_c31628_571, Tdurum_contig43995_611, Excalibur_c31707_302 demonstrated the highest effect hastening the maturity by about 9 days (Table 1).

Association of NAM-A1, Vrn-B3 and Vrn-B1 alleles with maturity time

In the course of a grain protein content study performed earlier, we had genotyped the NAM-A1 gene in the studied panel of wheat cultivars (Leonova et al., 2022). In this study, we compared the MT estimated for 11 years and NAM-A1 alleles. The results indicated that the NAM-A1 allele variation had no significant impact of on the MT (Fig. S1).

Allelic variation at the Vrn-B3, Ppd-D1, Vrn-A1, Vrn-B1 and Vrn-D1 loci was determined by Berezhnaya et al. (2021). The authors observed almost no variation at the Vrn-A1, Vrn-D1, and Ppd-D1 loci. Since Vrn-B1 and Vrn-B3 are polymorphic, their MT association was assessed. The results de/monstrated that the Vrn-B1 and Vrn-B3 alleles or their combination had no significant effect on the MT (Fig. S2).

Detection and prioritization of candidate genes

Significant SNPs spanned 21.9 Mb region of 7BL that included 216 genes. We estimated the functional annotation of these genes and their expression profiles in different tissues and developmental stages to select mostly appropriate candidate genes.

When evaluating the candidates’ expression patterns, we relied on information about the expression of known MT-associated genes. NAM1 are the key genes determining MT duration. The wild-type Nam-B1 allele has had a negative effect on yield traits, and it is rare in wheat cultivars (Hagenblad et al., 2012). Therefore, we couldn’t assess its expression pattern. The gene’s homoeologues, Nam-A1, and Nam-D1, are predominantly expressed in the flag leaf after the beginning of flowering, as reported on the website http://www.wheat-expression.com/.

To select the genes whose expression was specific to the flag leaf, the data obtained by Ramírez-González et al. (2018) were analyzed. The expression patterns of 216 genes in different tissues were assessed and the genes with preferable expression patterns were selected (File S1). There were twelve genes expressed in the flag leaf: TraesCS7B02G391800, TraesCS7B02G395600, TraesCS7B02G396600, TraesCS7B02G398700, TraesCS7B02G399700, TraesCS7B02G400500, TraesCS7B02G400600, TraesCS7B02G401600, TraesCS7B02G404300, TraesCS7B02G411800, TraesCS7B02G412200, TraesCS7B02G413400 (Table S5).

Expression estimation of candidate genes

To select the best reference genes for the plant developmental stages and conditions, the expression of seven previously published reference genes was evaluated. Analysis of the expression stability with the geNorm algorithm demonstrated that the best reference genes for wheat in anthesis phase were the genes encoding scaffold-associated regions (SAR) binding protein and methionine aminopeptidase 1 (MetAP 1) (Fig. S3).

To estimate the candidate gene expression, the most early- and late-maturing varieties having respective haplotypes CCGACTGACT and TAAGTAGTC within the QTL on 7B were analyzed. It was found that eight candidate genes TraesCS7B02G411800, TraesCS7B02G401600, TraesCS7B02G400600, TraesCS7B02G399700 TraesCS7B02G404300, TraesCS7B02G400500, TraesCS7B02G398700, TraesCS7B02G413400 were not expressed. Expression of three genes (TraesCS7B02G391800, TraesCS7B02G395600, TraesCS7B02G396600) was higher in early-maturing cultivars, but no significant difference in their expression level was found. On the contrary, the expression of TraesCS7B02G412200 was significantly lower in the early samples (p < 0.01) (Fig. 3).

Validation of association analysis results

Evaluation of the mapped population (Obskaya 2 × Novosibirskaya 15) showed that the offspring MT varied from 74.5 to 91.2 days. The SNP markers mapped on chromosomes 2A, 4A, 5B, and 7B (Table 1) were polymorphic between the parental varieties and the F_2:3 progeny, while the markers located on chromosomes 3B and 7A were found to be non-polymorphic according to the genotyping performed with an Illumina 15 K Wheat array. There was a significant difference in MT for wsnp_Ku_rep_c68565_67614479 and Excalibur_c53111_215 localized on chromosomes 4A and 5B, respectively (Fig. 4). No significant MT difference was shown found for SNP markers Tdurum_contig47508_250 and BS00023075_51 (chromosome 2A).

The markers on chromosome 7B were divided into two groups according to their genetic linkage and physical map location. Group I included four markers (Tdurum_contig43995_370, Kukri_c31628_571, Tdurum_contig43995_611, Tdurum_contig44993_359) and two haplotypes (TAAG and CCGA). Group II included seven markers (RAC875_c89312_61, IAAV8521, Tdurum_contig5360_329, IAAV6659, GENE-4442_121, Excalibur_c31707_302, BS00047083_51) and two haplotypes (TGAGTCA and CTGACTG). Comparison of the F_2:3 progeny showed that the TAAG and TGAGTCA genotypes matured later compared to the CCGA and CTGACTG ones (Fig. 4). These MT-related differences had a pronounced trend, despite the lack of significance. It should also be noted that in all cases, the alleles associated with later maturation were inherited from cv. Obskaya 2.

Association of the known genes with maturity time

To date, NAM1 are the only genes whose influence on the time of wheat ripening has been confirmed. It is known that the wild-type allele of the NAM-B1 gene is associated with earlier maturity, while it is almost not found in modern cultivars (Hagenblad et al., 2012). Also, the NAM-A1a and NAM-A1b alleles are specific for early and medium-early varieties (Alhabbar et al., 2018).

In previous study, analysis of the spring common wheat panel identified the NAM-A1a, NAM-A1c, and NAM-A1d alleles (Leonova et al., 2022). Since in our sample only five varieties carried the NAM-A1a allele, which together with NAM-A1b, is associated with early maturity, it is no wonder that no significant relationship between the allelic state of the NAM-A1 gene and MT was detected.

Vernalisation and Photoperiod genes are key factors of wheat development and some works have demonstrated that they may determine not only heading and flowering times but maturity time as well (Shcherban, Emtseva & Efremova, 2012; Kamran et al., 2013; Zaitseva & Lemesh, 2015; Whittal et al., 2018).

The wheat samples in our study were analyzed to identify the alleles of the Vernalisation genes. The plants had different alleles of Vrn-B1 and Vrn-B3 genes, while Vrn-A1, Vrn-D1, and Ppd-D1 were not found to be polymorphic. However, we did not observe any significant impact of Vrn-B1 and Vrn-B3 on the heading time of the examined wheat plants.

Among the VRN-1 gene alleles, accelerated development was predominantly associated with the Vrn-A1a allele, which was present in the majority of the sampled varieties. Similarly, the dominant alleles of Vrn-B1 (Vrn-B1a and Vrn-B1c) detected in our panel of cultivars were unlikely to have a different effect on development.

In our investigation of the Vrn-B3 gene, we identified two alleles: vrn-B3, which is the recessive variant typically found in late cultivars, and the Vrn-B3e allele, which delays heading time by 1.5 days when compared to vrn-B3.

The ANOVA results didn’t demonstrate a significant impact of the Vrn-B3 alleles on the maturity time. We also determined that under the conditions of the experiment, the main contribution to maturity time inputed the period from heading to maturity. Based on this finding, it can be suggested that Vrn-B3e does not appear to alter MT (maturity time) in spring wheat cultivars.

QTL colocalization with previous studies

In this study significant MT-associated SNPs were found on chromosomes 2A, 3B, 4A, 7A, and 7B. The most significant markers were located on the 7B chromosome. Previously, Kulwal et al. (2003) had reported about a QTL on the long arm of the 7B chromosome to be associated with MT. The QTL was tightly linked with RFLP marker XksuD2 located on 7BL (Paillard et al., 2003). The Blast analysis of the XksuD2 primers specified that it was located in the 683,517,074–683,517,093 region of 7BL, a position to be very close to the MTAs detected on 7BL in our study.

The second locus coincided with the previously detected region was positioned on the short arm of the 7A chromosome. Its MTA was associated with SNP marker Tdurum_contig20214_279 positioned in 61,816,447. The location of QMat.dms-7A.2 (Perez-Lara et al., 2016) was associated with RAC875_c14982_577 (position 59,045,900). So, these two loci may be assumed to be very close to each other and should be the same. However, MT-associated loci had been previously detected on the same chromosomes: 2A (Wang et al., 2009), 3B (McCartney et al., 2005), 4A (McCartney et al., 2005; Chen et al., 2015; Perez-Lara et al., 2016), 7A (Zou et al., 2017). Comparative analysis of their location demonstrated that they were different.

Candidate genes

NAM1, a maturity time gene, is expressed after flowering begins in the flag leaf. Moreover, a number of studies have demonstrated the flag leaf to be main source of assimilates during the grain-filling period, which concurs with senescence onset (Ashraf & Bashir, 2003; Viljevac Vuletić & Španić, 2020). This is the reason why we focused on the candidate genes in the QTL that were most strongly expressed in the flag leaf.

Annotations of these twelve genes are summarized in Table S5. Most of them are predominantly expressed in flag leaf after the heading or flowering begins.

One of the genes showing strong expression in the flag leaf at the flowering stage was TraesCS7B02G391800 (Fig. 5A). Three the most significant SNPs (Tdurum_contig43995_370, Kukri_c31628_571, Tdurum_contig43995_611) were linked with TraesCS7B02G391800 (Fig. 5B). These SNPs are either located in the last exon (Tdurum_contig43995_611) or in 3′UTR region (Tdurum_contig43995_370, Kukri_c31628_571). The gene encodes a BZIP domain-containing protein. BlastP for Arabidopsis sequences demonstrated that the protein encoded by TraesCS7B02G391800 was related to such bZIP transcription factors as bZIP9, BZO2H1 (bZIP10), bZIP25, bZIP63, and the bZIP9 protein known to regulate bZIP53, which, upon dimerization with bZIP10 or bZIP25 (both C-group bZIPs), may enhance maturity-gene expression (Alonso et al., 2009; Jain, Shah & Rishi, 2018). TraesCS7B02G391800 is ortologous to rice OsbZIP52 (Os06g0662200) that encodes a probable transcription factor that binds to the DNA specific sequence 5′-TGAGTCA-3′ found in seed storage protein gene promoters (PubMed:11133985).

However, qPCR revealed no significant difference of TraesCS7B02G391800 expression in the early and late maturing varieties, we suppose this gene to be a highly probable candidate.

TraesCS7B02G401600 is expressed in the flag leaf while anthesis and then in the embryo proper while the ripening stage. The gene is known as BRITTLE1 (BT1) and provides unidirectional transmembrane transport of ADP-glucose that is crucial for starch synthesis in cereal grain (Wang et al., 2019).

KnetMiner demonstrated that genes TraesCS7B02G391800 (bZIP9), TraesCS7B02G401600 (BT1/BT2) and TraesCS7B02G413400 may interact each other (Fig. 5C), which allows us to assume that in the identified locus, several genes might affect MT thanks to their interaction with each other.

TraesCS7B02G412200 (encode photosystem II reaction center) is expressed in the flag leaf after the flag-leaf stage begins. It was the only gene that demonstrated significant expression difference between early- and late-maturing variety groups in our qPCR experiment. Previous studies had demonstrated that the genes for chlorophyll binding, chloroplast, chloroplast organization, photosystem I, and photosystem II were down-regulated in early senescing varieties (Yoshida, 2003; Zhang et al., 2018).

KnetMiner analysis of the gene networks involving TraesCS7B02G412200 revealed its relation to CGA1 (CYTOKININ-RESPONSIVE GATA TRANSCRIPTION FACTOR) (Fig. S4). Cytokinins are key regulators of plant growth and development, including senescence (Argueso, Raines & Kieber, 2010). And CGA1 overexpression is known to be associated with delayed senescence and reduced grain filling (Hudson et al., 2013). The decay of chloroplasts produces most of the nitrogen reassimilated to the seeds (Masclaux-Daubresse, Reisdorf-Cren & Orsel, 2008). Hudson et al. (2013) have hypothesized that expression of CGA1 results in increased activity of chloroplasts, that in turn causes the delayed senescence, inhibits nutrient remobilization into seeds and prevents proper grain filling, which makes it possible to assume that TraesCS7B02G412200 and the genes associated with it can indirectly influence the maturity time through the regulation of photosynthesis and chloroplast functioning.