The study was conducted using a diverse collection of 225 Oryza sativa accessions, consisting of 83 accessions from the Mini Core Collection of Huazhong Agricultural University and 142 accessions from the 3K Rice Genomes Project. These accessions originated from various parts of the world and encompassed different subpopulations, which can complete their reproductive cycle in Wuhan. Table S1 provides the details of the accessions, including their names and countries of origin.

The experiment was conducted in the field of Huazhong Agricultural University, Wuhan, China. About 200 g of seeds were sown on May 15^th of 2018 and May 18^th of 2019. 30-day-old seedlings were transplanted into 1 m × 2 m plots, with one plant per hill at a spacing of 0.20 × 0.25 m. Fertilizers applied to all plots were 180 kg N ha^-1, 60 kg P₂O₅ ha^-1, and 120 kg K₂O ha^-1. The plots received standard management practices, including irrigation, fertilization, and disease and pest control. Figure S1 shows the weather data for the whole growing season.

To unify the data from two years, we used the lmer function within the lme4 package. The phenotype data were modeled with a linear mixed model, where accession was the fixed effect and year and replication were the random effects, to calculate the BLUE (best linear unbiased estimator, fixed factor) values to be used in the GWAS analysis. The following formula was used to calculate the heritability:

, where, V G , V G E , V ϵ , r , and e represent the genetic variance, interaction variance between genotypes and environments, error variance, number of replicates within each environment, and number of environments, respectively. Data entry was done using MS Office, while analysis and processing were carried out using the R software (https://cran.r-project.org/).

Five plants of each accession in the middle of the plot were selected to investigate the chlorophyll content and fluorescence characteristics at the heading stage. Chlorophyll content and fluorescence parameters were measured in the middle part (1/3~2/3) of the flag leaves between 8:30 and 11:30 a.m. on a sunny day, using a portable chlorophyll fluorometer (MultispeQ v1.0) to obtain more reliable data in the field setting. The instrument was used with the protocol “Leaf Photosynthesis MultispeQ V1.0 no open/close” provided at https://www.photosynq.org/protocols/leaf-photosynthesis-multispeq-v1-0-no-open-close, which is a classic and by far the most utilized PhotosynQ Protocol for measuring many photosynthesis-related parameters in a short period of time. Due to insufficient dark adaptation during the measurement, our F v / F m parameter is not rigorous and can only reflect the plant’s state at the time of measurement. The SPAD of the leaf was calculated by measuring the absorbance at 650 nm and 940 nm, and k is the calibration coefficient obtained using MultispeQ calibration cards.

The eight fluorescence parameters were calculated based on the minimum fluorescence ( F o ), maximum fluorescence ( F m ), steady-state fluorescence ( F s ), F o ' and F m ' (the same as above but measured under light conditions), and photosynthetically active radiation (PAR). The calculation of Φ I I , Φ N O , Φ N P Q , and F v / F m parameters related to photosynthetic efficiency was carried out as follows:

The linear electron flow (LEF) was calculated as follows:

The q L and q P , which reflect the “Lake” model and “Puddle” model in Photosystem II Redox State, was calculated as follows:

N P Q t , an efficient parameter that reflects NPQ, was calculated without the need for complete relaxation of the quenching process. The calculation for N P Q t was as follows:

DNA was extracted from fresh leaves of field-grown plants using a modified CTAB method (Yan et al., 2008). Whole-genome DNA sequencing was performed on the Illumina HiSeq-2000 platform by Personalbio (Shanghai, China). (Andrews, 2010) (V0.11.9) was used for quality control of sequencing data, and paired-end 150 bp reads were mapped to the Nipponbare reference genome (https://www.ebi.ac.uk/ena/data/view/GCA_001433935.1) using BWA (V0.7.17) with the default parameters. After alignment, the genomic data were sorted using SAMtools (V1.9) and the sequencing reads were de-duplicated using SAMBAMBA (V0.8.2). Genomic variants (in GVCF format for each accession) were identified using the Genome Analysis Toolkit (GATK V4.3.0) software, with the HaplotypeCaller module and GVCF model. The raw variant sites were further filtered by Plink (V1.9), with genotype quality for each individual ≥ 10%. After genotype imputation using Beagle (V4.1), the minor allele frequency (MAF) was controlled to be ≥ 5%. The identified SNPs were further annotated using the ANNOVAR software (version 16-Jul-2017).

All 442,634 identified SNPs were used to build a phylogenetic tree and perform principal component analysis (PCA). The individual-based neighbor-joining (NJ) tree was constructed using the phylip (V3.697) and EvolView (http://www.evolgenius.info/), based on the p-distance and with 1,000 bootstrap replicates. PCA was conducted using the Plink (V1.9) with the command “–pca 10” to output the top 10 PCA results. Since the first three principal components are more representative, we utilized the top three PCA results in the subsequent GWAS analysis. To estimate the LD in our rice population, the squared correlation coefficient (r² ) between pairwise SNPs was computed using PopLDdecay (Zhang et al., 2019). The r² value was calculated for pairwise markers in a 1000-kb window and averaged across the whole genome. The “–cv” command of Admixture (V1.3.0) was used to calculate the cross-validation error for K = 2, 3, 4, and 5.

GWAS was performed using a mixed linear model (MLM) in the GEMMA (V0.98.1) package (Zoubarev et al., 2012). The matrix of pairwise genetic distances calculated by GEMMA was used as the variance-covariance matrix of random effects. The kinship matrix kin.sXX.txt was calculated using the command “-gk 2 -p Phenotype” and GWAS analysis was conducted using the command “-k kin.sXX.txt -lmm 1 -p Phenotype -c PCA”. Significant p-value thresholds P< 1.13 × 10^-7 (0.05/442,634) were set to control the genome-wide type 1 error rate, which was calculated by 0.05/n (total SNPs). PVE of 100 kb was filtered out before and after the peak signal. The Manhattan and quantile-quantile (QQ) plots of GWAS results were generated in R software (https://cran.r-project.org/).

To consistently estimate the genetic effect of Φ I I and SPAD, the genetic variants were selected according to the three assumptions in MR analysis, (i) the genetic variants were obtained from the results of GWAS associated with the single component trait at a genome-wide significant level (P < 1.13 × 10^-7); (ii) the genetic variants were not associated with any confounders; (iii) the genetic variants only affected SPAD through the Φ I I trait, not through other component traits (P > 0.05).

The MR Egger, Weighted Median, Inverse Variance Weighted, Simple Mode and Weighted Mode methods were used for MR analysis to assess the effect of Φ I I on SPAD, by summarizing the effects of multiple independent SNPs. In sensitivity analysis, the MR Egger method and Inverse Variance Weighted method were used for MR analysis. According to the results, leave-one-out analysis was supplemented. MR analysis was performed in R package TwoSampleMR (https://mrcieu.github.io/TwoSampleMR/).

Total RNA was separately extracted from each sample using an RN38 EASYspin plus Plant RNA kit (Aidlab Biotech, Beijing, China). RNA integrity was determined through the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States). The libraries were sequenced by Personalbio (Shanghai, China) with an Illumina HiSeq (Illumina, CA, United States) system. To ensure the accuracy, reads with more than 10% N bases and low-quality reads with Q ≤ 20 and over 50% bases were excluded (Chen et al., 2018). The resulting clean reads were mapped to the Nipponbare reference genome using Tophat2 (Kim et al., 2013). Gene expression was then calculated by counting the number of mapped clean reads for each gene normalized into Fragments Per Kilobase of transcript sequence per Millions (FPKM).

DESeq2 R package (Love et al., 2014) was used for multiple testing correction of DEGs, and the false discovery rate (FDR) was calculated through the Benjamini and Hochberg’s method. DEGs were defined as genes exhibiting at least a 2-fold difference in expression, with |log2FoldChange(L2FC)| > 1, and P < 0.05. The Pheatmap R package (https://www.rdocumentation.org/packages/pheatmap/) performs bidirectional clustering analysis on the union of all DEGs and samples in all comparison groups. The clusterProfiler R package (Yu et al., 2012) was used to perform GO enrichment analysis for DEGs, with the p-value adjusted through the Benjamini and Hochberg’s method and a P < 0.05 selected as the threshold for determining significant GO terms. For all samples, PCA was carried out to explain their interrelationship. Blast2GO (Conesa et al., 2005) was used for DEG annotation and functional prediction.

The filtering generated a total of 632.17 GB of high-quality reads, which were mapped to the Nipponbare reference genome, with an average success rate of 94.4% ( Table S2 ) and an average sequencing depth of 16.5-fold ( Table S2 ). A total of 9,989,556 SNPs were identified on 12 chromosomes from the mapping, with the highest and lowest density of SNPs being detected on chromosome 11 and chromosome 03, respectively, and the average marker density was 27.10 SNPs/kb ( Table S3 ). After filtering out SNPs with a low genotyping rate using PLINK, gene imputation was performed using Beagle. A final set of 442,634 SNP markers with a MAF greater than 0.05 was retained for GWAS analysis ( Table S4 ).

PCA ( Figure 1B ) divided the population into four groups, which is consistent with the results of phylogenetic tree ( Figure 1A ), and the K value was considered as the number of subgroups with the lowest cross-validation (CV) error ( Figures 1D , S2 ). Figure 1C shows the average linkage disequilibrium (LD) decay in the whole genome. Figure 1 indicated that these materials could be divided into four groups with some genetic differences from each other.

In order to reveal the fluorescence characteristics of 225 rice accessions, we evaluated the chlorophyll content (SPAD) and eight fluorescence parameters ( N P Q t , Φ I I , Φ N O , Φ N P Q , LEF, F v / F m , q L and q P ) in two years (2018 and 2019). Figure 2A shows the distribution of each trait. We used the BLUE value to combine the results of the two years, and performed descriptive statistical analysis ( Table S10 ). The box plots showed differences (Wilcox test for one group) in fluorescence parameters among the four groups classified by PCA ( Figure 2B ), where N P Q t (P = 0.025), Φ N P Q (P = 0.035), SPAD (P = 6.7 × 10^-9), and F v / F m (P = 0.027) exhibited significant differences among different PCA groups, indicating that the phenotypic differences in SPAD, N P Q t , Φ N P Q , and F v / F m among different accessions were affected by the genetic background.

We calculated the Pearson correlation coefficients to investigate the relationship between different fluorescence characteristics. As expected, SPAD, an indicator of chlorophyll content, was negatively correlated with Φ I I , q L , and q P (R = -0.21, -0.20, -0.21; P< 0.05, respectively). N P Q t , which can reflect non-photochemical quenching, showed a significant positive correlation with Φ N P Q (R = 0.90; P< 0.05). Φ I I , an indicator of photochemical efficiency, exhibited significant positive correlations with F v / F m , q L and q P (R = 0.41, 0.94, 0.99; P<0.05, respectively). Additionally, Φ N O showed a significant positive correlation with F v / F m , but negative correlations with Φ N P Q , LEF, q L and q P . LEF exhibited significant positive correlations with N P Q t and Φ N P Q , but negative correlations with Φ I I , Φ N O , and F v / F m ( Figure 3A ).

The heritability of traits is a key parameter in breeding selection (Nirmaladevi et al., 2015; Roy and Shil, 2020). Here, SPAD, N P Q t , Φ I I , Φ N O , Φ N P Q , LEF, F v / F m , q L , and q P exhibited different patterns of heritability, ranging from 0.06 to 0.97 ( Figure 3B ). The heritability of SPAD, N P Q t , and LEF was greater than 0.4, while that of Φ N P Q , F v / F m , Φ N O , Φ I I , and q L was below 0.4. SPAD had the highest heritability of 0.97. These results indicated that SPAD, N P Q t , and LEF are greatly influenced by genetic factors.

We conducted a GWAS using the MLM method implemented in GEMMA software and analyzed the final set of 442,634 SNPs. Q and K, which can represent the population structure and kinship, were included in the MLM model to prevent spurious associations, with a significance threshold of P< 1.13 × 10^-7. By integrating the Manhattan plots for rice chlorophyll fluorescence traits ( Figures 4A–L ) and LD decay rates of 12 chromosomes in 225 rice accessions ( Figure 1C ), and based on the LD coefficient decreasing to half of its maximum at a distance of 1 kb, we selected target intervals at 2 kb upstream and downstream of the SNP, and finally identified 31 significantly associated loci. These loci included 78 SNPs associated with SPAD, N P Q t , Φ I I , Φ N P Q , q L and q P , which comprised 3, 64, 8, 18, 26, and 5 SNPs, respectively ( Table S5 ). Moreover, clear co-localization was observed between Φ I I and Φ N P Q and between q L and q P ( Figure S3 ), and the co-localization results of Φ I I , q L , and q P were annotated in the Manhattan plot ( Figures 4C, I, K ). Based on functional analysis of genes in LD regions, a total of 17 candidate genes were identified for chlorophyll fluorescence characteristics ( Table S6 ). No significant SNPs were found for  Φ N O , LEF, and F v / F m ( Figure S4 ).

In section 3.2, we observed a negative correlation between Φ I I and SPAD. To comply with the requirements of MR analysis, we included 65 Φ I I loci that reached genome-wide significance (P< 1.13 × 10^-7) in the GWAS analysis. These loci, which exhibited negative genetic effects on SPAD, were consistently observed across five analytical methods ( Figure 5 ): MR Egger (Beta = -40.47; P< 0.05), Weighted Median (Beta = -20.43; P< 0.05), Inverse Variance Weighted (Beta = -19.94; P< 0.05), Simple Mode (Beta = -30.35; P< 0.05), and Weighted Mode (Beta = -29.69; P< 0.05) ( Table S8 ). In sensitivity analysis, homogeneity statistics showed that the effect sizes of the studied loci were homogeneous in MR Egger (P< 0.05) and Inverse Variance Weighted (P< 0.05) methods ( Table S8 ). As the Horizontal pleiotropy analysis result was insignificant (intercept = 2.10; P > 0.05) ( Table S8 ), we conducted a leave-one-out analysis on the 65 SNPs ( Table S9 ). The results further confirmed the negative effect of Φ I I on SPAD.

To further investigate the genetic basis for chlorophyll fluorescence characteristics in rice, two rice varieties with significant differences in Φ I I and SPAD, namely D062 (High Φ I I type, H) and D133 (Low Φ I I type, L), were selected from the population for further analysis. The Φ I I and SPAD of H and L are presented in Figures 6A–C . We collected flag leaves at the heading stage for RNA sequencing with three biological replicates for each accession. Finally, a total of 263,471,192 reads with a Q30 score of 92.33% were generated. Among these reads, 3.70–4.38% were multiply mapped, while 95.62–96.30% were uniquely mapped to the reference genome ( Table S7 ). To evaluate the data reliability, correlation analysis ( Figure 6D ) and cluster analysis ( Figure 6E ) were conducted. DEGs between the two varieties (H and L) were identified, including 1,434 upregulated genes and 932 downregulated genes ( Figure 6F ).

Expression clustering can identify the unknown biological connections between genes. Both H and L type had good correlations within the group, indicating that DEGs in different groups may have specific connections with certain biological processes, metabolisms, and signaling pathways ( Figure 7A ). To uncover the functions of 2,366 DEGs, Gene Ontology (GO) enrichment analysis was conducted, and the DEGs were classified based on their molecular function (MF), biological process (BP), and cellular component (CC). The top five GO terms with the smallest p-values, namely the most significant enrichments, were selected and presented for each category. For the MF category, the top five enriched GO terms were protein phosphorylation (GO:0006468), phosphorylation (GO:0016310), adenyl ribonucleotide binding (GO:0032559), adenyl nucleotide binding (GO:0030554), and protein kinase activity (GO:0004672). Based on the analysis results, these DEGs were likely involved in a host of biochemical reactions necessary for kinase activity, nucleotide binding, and phosphorylation ( Figure 6C ). For the BP category, the top five enriched GO terms were cell surface receptor signaling pathway (GO:0007166), response to stimulus (GO:0050896), protein phosphorylation (GO:0006468), phosphate-containing compound metabolic process (GO:0006796), and phosphorylation (GO:0016310). The results suggested that these DEGs may have crucial functions in cellular signaling, response to environmental stimuli, and metabolic processes ( Figure 6C ). For the CC category, the top five enriched GO terms were intrinsic component of membrane (GO:0031224), integral component of membrane (GO:0016021), plasma membrane (GO:0005886), membrane part (GO:0044425), and membrane (GO:0016020), indicating that these DEGs may be involved in various cellular membrane-related functions ( Figure 7B ). Further scrutiny of the DEGs indicated their potential involvement in regulating the photosynthetic performance of rice, including N P Q t , Φ I I , and SPAD. Therefore, these genes represent valuable research targets for further investigation and potential avenues for crop improvement.

To further confirm the candidate genes, we validated the genes related to rice chlorophyll fluorescence characteristics by combining GWAS significant regions, LD decay, DEGs, and gene annotation. The Os03g0583000 and Os06g0587200 genes were found to be located on the SNPs identified in the GWAS, and showed significant differences in expression levels (|L2FC| > 1, P< 0.05) between the two rice varieties. Therefore, these two genes were considered as the most likely candidate genes. Blast2GO annotations revealed that Os03g0583000 was a peroxisomal protein, and Os06g0587200 was a protein kinase containing a catalytic domain ( Table 1 ).

We plotted the genetic structure of the two candidate genes ( Figure 8G ), and haplotype analysis showed that six SNPs in the promoter of Os03g0583000 formed two haplotypes ( Figure 8G ). The inbred lines carrying haplotype 1 had significantly lower Φ I I and q P values while significantly higher SPAD values than those carrying haplotype 2 ( Figures 8A–C ). In addition, transcriptome analysis showed that the H type had a significantly lower FPKM value of Os03g0583000 than the L type ( Figure 8E ), indicating that Os03g0583000 was the most likely candidate gene for Φ I I and q P . The genetic variations at the identified SNP loci were also found to affect SPAD. In addition, three SNPs in the exon region of Os06g0587200 formed two haplotypes ( Figure 8G ). The inbred lines carrying haplotype 1 had significantly lower N P Q t than those carrying haplotype 2 ( Figure 8D ). Transcriptome analysis showed that the H type had a significantly higher FPKM value of Os06g0587200 than the L type ( Figure 8F ), and the H type belonged to haplotype 2, whereas the L type belonged to haplotype 1, indicating that Os06g0587200 was the most likely candidate gene for N P Q t .