Genome-wide association study reveals genetic basis and candidate genes for chlorophyll content of leaves in maize (Zea mays L.)

Plant materials

The association panel utilized in this study included 334 maize lines provided by the Maize Institute of Sichuan Agricultural University. These lines were collected from the breeding program of Southwest China and consist of tropical, non-stiff stalk (NSS), stiff stalk (SS), and other unique germplasms (Zhang et al., 2016; Table S1). The majority of these accessions belong to the mid to late maturity group, with maturity periods ranging from 100 to 125 days. The panel was rigorously evaluated across three distinct environments with significant differences in climatic conditions, namely Chongzhou (CZ, Sichuan Province; 30.30° N, 103.07° E) and Ya’an (YA, Sichuan Province; 29.59° N, 102.57° E) in 2021, as well as Xishuangbanna (XSBN, Yunnan Province; 22.0° N, 100.79° E) in 2022. The evaluation employed a completely randomized design with three replicates for consistency and reliability. Each line was grown in a single row with row length of 3 m and row distance of 0.7 m. A standard corn management practices were applied during the cultivation.

Phenotypic data collection and analysis

The SPAD values were collected from five plants with consistent growth condition of each line to symbolize the chlorophyll contents using SPAD 502 Plus Chlorophyll Meter (a handheld SPAD instrument). At the seedling stage (30 days after sowing) and the grain filling stage (5 days after pollination), the CC was measured at the middle parts of the fifth leaf and the ear leaf, respectively. Each plant was measured three times, after which the mean value was recorded as the leaf CC. The descriptive statistics of CCFSS and CCEFS in each environment, including mean, maximum (Max), minimum (Min), standard deviation (SD), coefficient of variation (CV), skewness, and kurtosis of CC were analyzed using the psych package (Revelle, 2024) in R 4.4.2 (R Core Team, 2024). To evaluate multi-environment experimental data, the BLUP values were computed using a liner mixed model for the estimation of random effects with lme4 R package. The ANOVA analysis was performed to calculate the variance components of each trait, including genotypes (G), environments (E), and interactions between genotype and environment (G × E). The broad-sense heritability (H²) was estimated using the following formula (Knapp, 1986):

$$H^2={\displaystyle \frac{{\sigma }_g^2}{{\sigma }_g^2+{\displaystyle \frac{{\sigma }_{ge}^2}{e}+{\displaystyle \frac{{\sigma }_e^2}{er}}}}}$$

Here, ${\sigma }_g^2$, ${\sigma }_e^2$, and ${\sigma }_{ge}^2$ represent genetic variance, residual error variance, and the variance of genetic × environmental interaction, respectively. The letters “e” and “r” denote the number of environments and the number of independent replicates, respectively.

GWAS

The genotype of the association panel was genotyped using the Maize SNP50 BeadChip, which consisted of 56,110 SNPs (Zhang et al., 2016). Those SNPs with minor allele frequency (MAF) < 0.05, missing rate > 20%, or heterozygosity > 20% were considered as low-quality variations, which had been filtered out. Then, a total of 43,728 high-quality SNPs were used for GWAS based on MLM model in the GEMMA package. Meanwhile, population structure (Q = 6) reported in the previous study was used as a covariate in the GWAS model (Zhang et al., 2016). The SNPs were considered as trait-associated QTNs with the p-value less than 1/N = 4.49 × 10⁻⁵ (N = 22,277). The N denotes the effective marker number of independent tests, which was calculated using simpleM function in R 4.4.2 package. The quantile–quantile (Q-Q) and Manhattan plots for GWAS were generated by using the CMplot function in R 4.4.2 package (https://github.com/YinLiLin/CMplot). The PVE of each significant associated SNP was calculated according to the formula as follows (Liu et al., 2024):

$$PVE={\displaystyle \frac{2{\hat{\beta }}^{2}MAF\left(1-MAF\right)}{2{\hat{\beta }}^{2}MAF\left(1-MAF\right)+{\left(se\left(\hat{\beta }\right)\right)}^{2}2N\ast MAF\left(1-MAF\right)}}$$where $\hat{\beta }$ is the effect estimate of genetic variants, MAF is the minor allele frequency of genetic variants, N is the sample size, and $se\left(\hat{\beta }\right)$ is the standard error of the effect. The raw phenotype and genotype data for maize lines are available in the figshare database with the accession link of https://doi.org/10.6084/m9.figshare.26355523.v1.

Analysis of candidate genes

According to the LD decay of this panel, gene models located within the 220 kb flanking regions of all trait-associated QTNs were identified as potential trait-related candidate genes. The function descriptions, gene ontology (GO) terms, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of these genes were annotated based on the Annotation database of maize B73 RefGen v2 in MaizeGDB (https://www.maizegdb.org/, accessed 25 May 2024). Enrichment analysis of GO and KEGG were performed using an online tool OmicShare (https://www.omicshare.com/).

Gene-based association analysis

The variations within the gene bodies and 2,000 bp upstream regions of candidate genes were obtained by DNA re-sequencing from 77 maize lines (Liu et al., 2024). The MLM model was tested to detect key variations associated with CC of maize at the p-value of 0.05. The LD decay between pairwise SNPs was calculated using LDBlockShow software (version 1.4) (Dong et al., 2021).

Gene expression patterns

The expression levels of candidate genes of distinct tissues at different development stages in maize were obtained from a previous study (Stelpflug et al., 2016). The heatmap was drawn using the pheatmap function in R 4.4.2 package (Kolde, 2019; R Core Team, 2024).

Statistical analysis

The two-sided t-test of the CC-related traits between two types of haplotypes was performed in Excel 2021. Box plots were created using the R ggplot2 package (https://ggplot2.tidyverse.org/).

Phenotypic variation

Phenotypic data of 334 maize inbred lines were collected across three experiment environments (Chongzhou, CZ; Xishuangbanna, XSBN; and Ya’an, YA;). The values of mean, Max, Min, SD, CV, skewness, and kurtosis of CC showed significant variability (Table 1 and Fig. 1). The SPAD values of CCFSS varied between 29.78–58.54, 28.89–53.56, and 31.53–55.15 across CZ, XSBN, and YA environments, respectively, with the CVs of 10.10%, 10.25%, and 11.62%. While, the range SPAD values of CCEFS in CZ, XSBN, and YA were 42.70–66.09, 34.67–68.40, and 38.55–64.84 respectively, with the CVs of 7.72%, 7.45%, and 7.78%. In addition, the BLUP values of CCFSS and CCEFS were calculated to eliminate the environmental deviation, with an average of 41.54 (ranging from 34.93 to 50.39) and 53.69 (ranging from 45.58 to 62.13), separately. Moreover, the absolute values of skewness and kurtosis for CCFSS and CCEFS across all environments and BLUP were less than 1.0 except for kurtosis of CCEFS in XSBN, indicating that CC followed normal distributions, and was controlled by numerous genes. Besides, the SPAD values shown significant positive correlations between any two environments or BLUP, indicating that there is a certain connection between the chlorophyll content of different development stages (Fig. S1). The values of broad-sense heritability (H²) for CCFSS and CCEFS were 70.84% and 78.99%, respectively, which confirmed that the CC of leaves were mainly controlled by genetic factors (Table 1).

QTNs associated with CC by GWAS

We used a MLM method with a threshold p-value of 4.49 × 10⁻⁵ to identify CC-related genetic loci. In total, 15 CC-related SNPs were detected. Among them, 11 SNPs were associated with CCFSS, and four were associated with CCEFS (Figs. 2, 3 and Table S2). The PVE values for these SNPs were 4.98–7.59%, indicating that chlorophyll content in maize was controlled by multiple mini-effect genetic loci. Notably, four SNPs were identified as co-located loci, which were detected in at least two different environments (including BLUP). These common loci indicate that they have more stable genetic effects, which should be attentioned in further studies. Three of the co-located SNPs (PZE-101214133, PZE-106069023, and PZsE-108003930) linked to CCFSS were suited in chromosomes 1, 6, and 8, respectively. One SNP (SYN23593), associated with CCEFS, was located in chromosome 7. For PZE-101214133, PZE-106069023, and PZE-108003930, the mean value of CC of all germplasms with the minor allele was significantly (p-value < 0.01) higher than that of germplasms with the major allele across the population (Table 2). However, the opposite performance was observed in SYN23593, where the average CC value of germplasms with the major allele was significantly higher than that of germplasms with the minor allele. (Table 2). Specifically, for CCFSS, the most significantly associated marker was PZE-101222689 (p-value = 4.68 × 10⁻⁷), located on chromosome 1 and explaining 7.37% of the phenotypic variation (Table S2). This SNP was located in the first intron of the gene GRMZM2G416388, which encodes a cystathionine beta-synthase (CBS) family protein. CBS proteins were known as energy sensors that regulate protein activities via their adenosyl ligand binding capacity. Specifically, in Arabidopsis, CBSX2 has been shown to inhibit the activities of m-type thioredoxins (TRXs) toward two chloroplast TRX-related targets, thus regulating plant growth (Baudry et al., 2022). For CCEFS, the most significantly associated marker was SYN23593 (p-value = 2.95 × 10⁻⁷), with a PVE of 7.59 (Table S2). This marker was located in the 3′-UTR of the gene GRMZM2G057296, which encodes a pectin lyase-like superfamily protein (PEL). In rice, the PEL gene takes part in the regulation of plant growth and leaf senescence (Leng et al., 2017).

We then analyzed the superior allele ratios in 30 maize lines widely utilized in the Southwest of China to evaluate the superior allele application of the 15 significant SNPs during maize breeding. Since higher CC is important for crop production, we considered the allele with the positive effect as the superior allele. Conversely, the alleles associated with lower SPAD values were designated as the inferior alleles. The superior allele percentage of all CC-related SNPs ranged from 0% (PZE-102047226, PZE-103162301, and PZE-104118967) to 70% (PZE-101222689) (Fig. 4). Notably, only three SNPs (PZE-101222689, PZE-106067897, and PUT-163a-149111517-965) had a superior allele ratio greater than 50%. However, eight SNPs (PZE-101068676, PZE-101214118, PZE-101214133, PZE-102047226, PZE-103162301, PZE-104118967, PZE-107016552, and PZE-108003930) possessed superior allele ratios of <10%, implying that these eight SNPs should be given priority consideration in molecular marker-assisted breeding to modify the CC of maize. Especially, two of the four co-located SNPs (PZE-101214133 and PZE-108003930) respectively existed in three and two maize lines, which should be paid more attentions. Among the 30 elite inbred lines, the maize line Mian723 owned the greatest number of superior alleles, possessing higher CCFSS and CCEFS of 47.79 and 62.13 in BLUP, respectively. The line F06 had only one superior allele, with the lower CCFSS and CCEFS of 38.03 and 49.98, respectively. Therefore, Mian723 and F06 had great utility values to display higher CC by providing or integrating more superior alleles.

Analysis of candidate genes

To further select the potential CC-related candidate genes, we searched the genes within 220-kb flanking regions of the 15 significant genetic loci identified by GWAS. In total, 177 candidate genes were detected, of which 131 genes had functional annotations (Table S3). Among these genes, 16 common genes were situated in the LD regions of SNPs PZE-101214118 and PZE-101214133, which were close to each other (<2,000 bp). Similarly, two genes simultaneously located in the LD region of PUT-163a-149111517-965 and PZE-106069023. According to the annotations, GRMZM5G820904 encodes a translocon at the outer envelope membrane of chloroplasts 75-III protein (TOC 75-III). As previously reported, TOC gene initiated the import process of thousands of nuclear precursor proteins, which are crucial for chloroplast formation, plant growth and development (Richardson et al., 2014). Another gene GRMZM2G164084 encodes an RNA polymerase sigma factor, which is essential to life and controls the process of transcription (Borukhov & Nudler, 2008). While the RNA polymerase sigma factor controls all transcription initiation steps and the stimulation of the primary steps in RNA synthesis (Vishwakarma & Brodolin, 2020). In Arabidopsis, gene SIG2 takes part in the transcription of several chloroplast tRNA genes possibly couples translation and pigment synthesis in chloroplasts (Kanamaru & Tanaka, 2004). The E1-E2 ATPase encoded by GRMZM2G324462 is also known as P-type ATPase. In Arabidopsis, a P-type ATPase involved in regulating the expression of a downstream gene ALA10, impacts the fatty acyl composition of chloroplast phosphatidylcholine, changing chlorophyll contents (Botella et al., 2016). GRMZM5G883222 encodes a phosphatidylinositol-4-phosphate 5-kinase family protein, which is involved in the initiation of chloroplast division by fusing a part of sequence of the prokaryotic FtsZ (a prokaryotic homolog of tubulin) (Shimada et al., 2004). In addition, the four co-located SNPs, PZE-101214133, PZE-106069023, PZE-108003930, and SYN23593 harbored 16, 9, 11, and 6 genes, respectively. Notably, the phosphate transporter encoded by GRMZM2G064657 (PZE-106069023), which was confirmed to affect the chlorophyll contents by regulating the phosphate acquisition in Cucumber (Naureen et al., 2018). The pentatricopeptide repeat (PPR) superfamily protein encoded by GRMZM2G071162 (PZE-106069023) is involved in the post-transcriptional regulation of chloroplast genes, and effects on the biogenesis and functioning of chloroplasts (Wang et al., 2021). These results further demonstrate that the candidate genes were potentially associated with CC.

To further reveal the function of these genes, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. As results, GO terms of “peroxidase activity”, “oxidoreductase activity”, and “inositol phosphate biosynthetic process” and pathways of “phenylpropanoid biosynthesis”, “MAPK signaling pathway”, and “inositol phosphate metabolism” were significantly enriched (p-value < 0.05), indicating that these genes hold significant potential for further research and exploration (Fig. 5, Tables S4, S5). Specially, there are 19 and 22 genes enriched in top 20 GO terms and KEGG pathways respectively, including eight genes that are previously identified as common located genes.

Gene-based association analysis revealed loci affecting maize CC

To further reveal the key variation loci affecting CC in maize, we conducted gene-based association analyses for the eight hub candidate genes using 77 randomly selected lines from the maize association panel. A total of 678 high-quality variations (542 SNPs, 69 insertions, and 67 deletions) located in the gene regions and their 2,000 bp upstream were obtained by DNA re-sequence (Table S6; Liu et al., 2024). Using the MLM model, a total of six variations (four SNPs, one insertion, and one deletion) from four genes (GRMZM2G037152, GRMZM5G816561, GRMZM2G324462, and GRMZM2G064657) were significantly (p-value < 0.01) associated with maize CC (Table S7). Among them, SNP-1-269340061 was situated in the first exon of gene GRMZM2G037152 and annotated as a missense variant, which probably alter the corresponding protein sequence. SNP-2-16709401 was annotated as a splice region variant, and located in the first intron on the gene GRMZM5G816561. SNP-6-126152668, DEL-6-126298795, SNP-6-126298598, and INS-6-126298833 were all located in the upstream region of the genes GRMZM2G324462 and GRMZM2G064657. GRMZM2G037152 encodes a GNS1/SUR4 membrane family protein, which was found to be associated with photosynthetic metabolism and biological processes of photomorphogenesis (Ma, Ding & Li, 2017). Cytochrome P450 encoded by GRMZM5G816561 is one of the most prominent families of oxidoreductases class enzymes. It catalyzes NADPH- and/or O2-mediated hydroxylation reactions in primary and secondary metabolism in various species (Xu, Wang & Guo, 2015; Chakraborty et al., 2023). GRMZM2G324462 and GRMZM2G064657 encodes an ATPase E1-E2 type family protein and a phosphate transporter protein, respectively. As mentioned above, the E1-E2 ATPases also known as P-type ATPases, which is related to the regulation of the contents of chlorophyll in Arabidopsis (Botella et al., 2016). The phosphate transporter could regulate phosphate homeostasis and photosynthesis in chlamydomonas (Tóth et al., 2024) and share similar features with chloroplast transit peptides in Arabidopsis (Versaw & Harrison, 2002). In particular, the expression level of this gene shown a significant difference between a yellow-green leaf mutant of maize and its wild type (Li et al., 2021). Together, these results further demonstrate that the candidate genes were potentially causal genes affecting maize CC.

According to these significant variants, we constructed haplotypes of each gene. For GRMZM2G037152, the corresponding SNP (SNP-1-269340061) with allele of A/C divided 77 maize lines into two groups. The mean values of both CCFSS and CCEFS traits in A-type lines were higher than those in C-type maize lines in all environments, except for CCEFS in YA and BLUP. Likewise, two alleles (A/C) of SNP-2-16709401 for GRMZM5G816561 and two alleles (A/C) of SNP-6-126152668 for GRMZM2G324462 both divided the 77 lines into two groups. The lines with Hap1 (A) for GRMZM5G816561 shown higher chlorophyl contents than that with Hap2 (C) in all environments. In contrast, lines with Hap1 (A) for GRMZM2G324462 shown lower chlorophyl contents than that with Hap2 (C) in all environments. In addition, based on three significant CC-associated variations (DEL-6-126298795, SNP-6-126298598, and INS-6-126298833) from GRMZM2G064657, two predominant haplotypes were classified (Hap1: AA-, Hap2: G-G). A t-test analysis revealed that the average CCFSS and CCEFS of Hap1 (AA-) lines were significantly (p-value < 0.05) higher than those of Hap2 (G-G) lines in all environments, except for CCFSS in YA (Fig. 6). Thus, the Hap1 (AA-) was regarded as a superior haplotype for improving chlorophyll contents in maize.

Expression patterns of the candidate genes

Based on a public reported gene atlas of maize tissues at different development stages (Stelpflug et al., 2016), we examined the expression profiles of all 177 candidate genes and constructed a heatmap (Fig. S2). The expression levels of these candidate genes were varied significantly in different tissues. Furthermore, we focused on eight hub CC-related genes, and found that GRMZM5G816561 showed a lower expression level in all tissues of different development stages except for some leaf-related tissues. Similarly, GRMZM2G064657 showed higher expression levels in leaves relative to other tissues. Universally, GRMZM5G810275, GRMZM2G467059, GRMZM2G073668, GRMZM2G002499, and GRMZM2G324462 had a relatively higher level of expression in whole development stages in all tissues (Fig. 7 and Table S8). Gene GRMZM2G037152 displays a relatively higher expression level in leaf-related and internode-related tissues. These results provide more information for revealing the mechanism of chlorophyll synthesis in further studies.