The association panel consisted of a diverse collection of 529 O. sativa accessions including both landraces and elite varieties. The details about the accessions, including accession name, country of origin, longitude and latitude origin, and subpopulation identity, have been reported previously (Wang et al., 2015; Xie et al., 2015) and are available at the RiceVarMap (http://ricevarmap.ncpgr.cn).

The sequencing, SNP identification, and imputation of the association panel have been reported in previous studies (Wang et al., 2015; Xie et al., 2015). The SNPs of the 529 O. sativa accessions are available at the RiceVarMap (http://ricevarmap.ncpgr.cn). There were 2,767,191, 1,857,866, 1,925,362, and 3,916,415 SNPs (minor allele frequency ≥0.05; the number of accessions with minor alleles ≥6) used in GWAS for groups of Indica, Japonica, Aus, and the whole panel, respectively. We performed GWAS using the linear mixed model (LMM) and the simple linear regression model (LR) provided by FaST-LMM program (Lippert et al., 2011). The population structure of Q matrix and kinship (K matrix) was taken into account as cofactor when performing association mapping using the LMM method. Using a method described by Li et al. (2012), the effective numbers of independent SNPs and suggestive thresholds were calculated. The effective numbers of independent SNPs were 571,843, 245,348, 235,880, and 757,578 for groups Indica, Japonica, Aus, and the whole panel, respectively. The suggestive P values set as thresholds to identify significant association signals by LMM were 1.75 × 10⁻⁶, 4.08 × 10⁻⁶, 4.24 × 10⁻⁶, and 1.32 × 10⁻⁶ for groups Indica, Japonica, Aus, and the whole panel, respectively. To obtain independent association signals, multiple SNPs exceeding the threshold in a 5 Mb region were clustered by r² of LD ≥0.25 and SNPs with the minimum P-value in a cluster were considered as lead SNPs.

Field trials were conducted in two years. The rice seeds were sown in the Experimental Station of Huazhong Agricultural University, Wuhan, China, in mid-May of 2012 and 2013. Seedlings about 25 days old were transplanted to the field. Each plot consisted of 4 rows with 10 plants each. The planting density was 16.5 cm between plants in a row, and the rows were 26 cm apart. Field management, including irrigation, fertilizer application and pest control, followed essentially the normal agricultural practice.

Three plants in the middle from the third row of each accession were subject to investigation for NPQ values before heading stage. Flag leaves were excised and floated on ion-exchanged water containing 0.01% Triton X-100 or 0.1% agar in centrifuge tubes. Leaves were dark-acclimated for at least 2 h. Chlorophyll fluorescence was measured with the PAM-2500 portable chlorophyll fluorometer (WALZ, Germany). Actinic lights were supplemented at a PPFD of 1,000 μmol m⁻² s⁻¹ for 5 min before measurements of Fm' to calculate parameter NPQ. The Stern-Volmer type nonphotochemical fluorescence quenching (Bilger and Björkman, 1990) was calculated by the following equation: NPQ = Fm/Fm'-1.

Three F2 populations were constructed for QTL mapping, with parents selected from the association panel. All of them were derived from inter-group crosses, an IndI by TrJ cross (Sadu-cho/Dourado Precoce), an IndI by TeJ cross (Zegu/Weiguo), and an Aus by TeJ cross (Dular/Akitakomachi). For each F2 population, 144 individuals were investigated for genotypic and phenotypic assay. Genotypic assay was performed using SSR markers, and measurement of NPQ values was conducted with the same method as described above. QTL analysis was conducted using composite interval mapping in QTL Cartographer (http://statgen.ncsu.edu/qtlcart/cartographer.html). The SSR markers used in the QTL analysis are presented in Table S4.

The genomic DNA extraction and sequencing were carried out according to the methods described by Du et al. (2011). The DNA sequences were assembled using SeqMan software and all polymorphisms were confirmed visually to ensure the accuracy. Sequences were inspected initially using Sequencher V5.4 program and edited manually using the SeqMan software. The genomic sequences were aligned by MAGE5. Nucleotide diversity and Tajima's D statistics were calculated using the DnaSP 5.0 program.

Total RNA was isolated from rice leaves using TRIzol Reagent (TransGen) according to the manufacturer's instructions. Approximately 3 μg of RNA sample was subjected to RNase-free DNaseI (Invitrogen) treatment and reverse transcribed using M-MLV Reverse Transcriptase (Invitrogen) with Oligo(dT)15. Quantitative RT–PCR was performed in a ViiA 7 Real-Time PCR system (Applied Biosystems) using FastStart Universal SYBR Green Master (Rox) superMIX (Roche). Ubiquitin was used as a reference gene in the qRT–PCR experiments. The measurements were obtained using the relative quantification method. The significant difference was analyzed statistically by student's t test or one-way analysis of variance (ANOVA). Primer pairs for qRT-PCR analysis are listed in Table S3.

To generate the RNAi constructs, the 425-bp fragment from the exon of OsPsbS1 was amplified by PCR using specific primers (Table S3), and cloned into the KpnI-BamHI and SpeI-SacI restriction sites of RNAi construct ds1301.

The CRISPR/Cas9 vector system for multiplex targeting of gene sites in monocot plants was kindly provided by Professor Yaoguang Liu (South China Agriculture University). In order to increase the efficiency of the targeting, we designed two target special sequences for each gene in exon regions. Briefly, the target site sequences were cloned into the sgRNA expression cassette which contains rice promoter, and then cloned into the BsaI restriction sites of vector pYLCRISPR/Cas9-MH. The primers used in constructing the sgRNA vectors for OsPsbS1 and OsPsbS2 are listed in Table S3.

The RNAi and CRISPR/Cas9 constructs were independently introduced into Agrobacterium tumefaciens strain EHA105, and then transformed into Zhonghua 11.

Histochemical staining of GUS activity in rice tissues was conducted essentially as described previously. Various tissues of transgenic-positive transformants (root, leaf, sheath, panicle, stem and mature seed) were incubated in GUS staining solution (50 mM sodium phosphate at pH 7.0, 10 mM Na₂-EDTA, 0.1% Triton X-100, 1 mg/ml X-Gluc, 100 μg/ml chloramphenicol, 1 mM potassium ferricyanide, 1 mM potassium ferrocyanide, and 20% methanol) at 37°C for 2–10 h after 15-min vacuum filtration. After GUS staining, the samples were incubated in 70% ethanol to remove chlorophyll and photographs were taken under a dissecting microscope (Leica MZFLIII).

As described by us previously (Wang et al., 2015; Xie et al., 2015), the whole association panel exhibited a distinct population structure, and was classified into 98 indica I (IndI), 105 indica II (IndII), 92 indica intermediate, 91 temperate japonica (TeJ), 44 tropical japonica (TrJ), 21 japonica intermediate, 46 Aus, and 14 VI (an intermediate group). The basic information of 529 accessions is available at the RiceVarMap (http://ricevarmap.ncpgr.cn).

The mean values and standard deviation of NPQ in different groups and all accessions in 2012 and 2013 (designated as NPQ_12 and NPQ_13, respectively) are presented in Table 1. A significant correlation was observed between NPQ_12 and NPQ_13 in the whole association population (r = 0.64). Large phenotypic variation was mainly detected between different subspecies in both years. Comparisons of the trait values between indica (consisting of IndI, IndII and indica intermediate) and japonica (consisting of TeJ, TrJ and japonica intermediate) revealed that japonica generally had significantly higher values in NPQ_12 (P = 2.79 × 10⁻⁴⁷) and NPQ_13 (P = 2.47 × 10⁻⁵⁹). Aus exhibited similar magnitude of NPQ values to indica (P > 0.05) in both two years, while much lower values than that in japonica (P = 7.68 × 10⁻¹¹ for NPQ_12 and P = 1.29 × 10⁻¹⁷ for NPQ_13). However, no statistically significant differences in NPQ values were observed between IndI and IndII or between TeJ and TrJ, with the only exception occurred between IndI and IndII in NPQ_13 (P = 0.00028).

Genome-wide association analyses were performed separately in the whole population and in the indica, japonica, and aus subpopulations for each year. A total of 33 significant association regions were identified at linear mixed models (LMM) (lead SNPs less than or around 300 kb were considered as caused by one common gene and counted as only one association region) for NPQ values. The details about these significant association signals are listed in Table 2. The quantile-quantile plots and Manhattan plots of LMM for NPQ_13 in the whole population are illustrated in Figure 1 as an example.

GWAS results over the two years were compared. Among them, two regions were commonly detected in both NPQ_12 and NPQ_13. Significant genetic heterogeneity was observed when comparing GWAS results performed in the different groups and in the whole population. Only three significant association regions were detected in at least two different groups (Table 2).

To search for candidate genes, the significant GWAS signals were compared with the positions of known or putative genes involved in NPQ regulation. We found that OsPsbS1 corresponds to one significant association region that was repeatedly detected and explained more than 40% of the variation in the whole population in both 2012 and 2013 (Table 2). However, no other previously identified NPQ key factors, such as genes related to xanthophyll cycle and lutein synthesis (Niyogi et al., 1998, 2001; Ruban, 2016), were found to be close to the detected significant association signals. It suggests the NPQ variation in rice may be controlled by a lot of new genes. Close to the GWAS peak signals, we found some candidate genes encoding light-induced protein 1-like, thylakoid lumenal protein, photosystem II (PSII) 10 kDa polypeptide, ATPase, ATP synthase gamma chain, and plastocyanin-like domain containing protein (Table 2). These candidate genes might be involved in processes of electron transport, buildup of proton gradient, PSII light harvesting antenna structure and light harvesting complex II (LHCII) rearrangements that trigger and regulate NPQ scenario (Ruban, 2016).

To investigate whether GWAS significant signals could be detected in bi-parental QTL mapping, we constructed three F2 mapping populations derived from inter-group crosses, with parents selected from the association panel. In the F2 population of IndI by TrJ (Sadu-cho/Dourado Precoce), one major QTL qSD1.1, explaining 29.4% of the phenotypic variation, located in the interval between RM212 and RM6504 on chromosome 1, was identified. In addition, a QTL qSD8.1, located in the interval between RM556 and RM80 on chromosome 8, also showed significant effects on NPQ values. The detailed information about these significant QTLs is presented in Table 3.

In the IndI by TeJ (Zegu/Weiguo) F2 population, similarly, one major QTL qZW1.1, explaining 33.6% of the phenotypic variation, located in the interval between RM8278 and RM315 on chromosome 1, was detected. Furthermore, a QTL qZW7.1, located in the interval between RM7571 and RM13.40 on chromosome 7, exhibited significant effects (Table 3).

In the F2 population of Aus by TeJ (Dular/Akitakomachi), one major QTL qDA1.1, explaining 37.0% of the phenotypic variation, located in the interval between RM212 and RM315 on chromosome 1, was also identified. Additionally, two QTLs qDA7.1 (located in the interval between RM6344 and RM420 on chromosome 7) and qDA9.1 (located in the interval between RM5657 and RM434 on chromosome 9), showed significant effects on NPQ values (Table 3).

After comparing the regions of these QTLs with GWAS signals, we found that all the major QTL detected on chromosome 1 corresponds to OsPsbS1, and other significant QTLs also co-localize with at least one GWAS signal (Table 3), indicating the GWAS results were quite reliable.

To characterize the nucleotide diversity of OsPsbS1, the whole genomic DNA sequences around OsPsbS1 from 480 accessions in the association panel were sequenced and analyzed. In total, 43 SNPs and 7 insertions and deletions (InDels) were detected in the aligned 6,997 basepairs, which include 4,022-bp promoter and 2975-bp coding regions (Table 4). Varied DNA polymorphisms were observed in different regions of the OsPsbS1 genome. We found no nucleotide polymorphic sites in exon-1, intron-1, and exon-2, and only one synonymous SNP in exon-3, indicating the structure and function of PsbS protein are highly conserved in rice. In the whole germplasm population, the pairwise nucleotide diversity parameter (π) and the level of the Watterson estimator (θ) in intron-2 and 3′UTR regions were much higher than that in the exons, promoter, and 5′UTR regions. Tajima's D values reached a significant positive level in the OsPsbS1 gene region (Table 4). Considering the strong population stratification, we also tested these parameters within the two subpopulations (indica and japonica group). The sequence diversities of indica in intron-2 and 3′UTR regions were also much higher than that in the exons, promoter, and 5′UTR regions. Compared with indica, lower levels of nucleotide diversity were observed in japonica. The Tajima's D values showed negative but not statistically significant values in both indica and japonica group (Table 4).

The sequencing results showed that the vast majority of japonica varieties has the 2,674-bp insertion fragment in the promoter region of OsPsbS1, compared with most indica and aus varieties. We further analyzed the effects of haplotypes based on 15 SNPs and 2 indels in the promoter region. For the strictness of statistical analysis, we only focused on those haplotypes shared by at least five accessions. The 480 accessions were accordingly divided into 7 haplotypes, among which hap1-3 contained the 2674-bp insertion, whereas hap4-7 didn't (Figure 2A). It was found that haplotypes with the 2674-bp insertion exhibited significantly higher values in NPQ_12 (Figure 2B) and NPQ_13 (Figure 2C) than other haplotypes without the insertion. We also compared the expression levels of OsPsbS1 in flag leaves of 49 randomly chosen accessions with the two major haplotypes, hap1 and hap4. Real-time PCR results showed that the expression level of OsPsbS1 in hap1 was much higher than that in hap4 (Figure 2D). The results suggest that variation of the OsPsbS1 expression levels might be the main cause of natural variation of NPQ capacity in the whole association panel.

To verify the function of OsPsbS1 in vivo, we first utilized RNA interference to silence the expression of OsPsbS1. The 425-bp fragment from the exon of OsPsbS1 was cloned into RNAi constructs ds1301, which were then transformed into Zhonghua 11, a japonica variety. Further real-time PCR analysis showed that the transcription of OsPsbS1 was drastically reduced in all transformants with much lower NPQ values (Figures 3A,B).

In the sequenced rice genome, OsPsbS2, another chlorophyll A-B binding protein gene, is found to be most similar to OsPsbS1 with a nucleotide identity of 91%. To distinguish the function of these two homologous genes, we then used the CRISPR-Cas9 system with constructs containing Cas9 and sgRNA and targeting the single gene OsPsbS1 or OsPsbS2 and both of the two genes. The constructs were transformed into a japonica variety Zhonghua 11. Two target special sequences were designed for each gene. The potential mutations were detected by sequencing the PCR amplification products containing the two targets regions. Indels mutations were confirmed in the target regions. Compared with the wild-type plants, the OsPsbS1 mutants and the double mutants exhibited much lower NPQ values (Figure 3C). However, no significant differences were found in NPQ values between the OsPsbS2 mutants and the wild-type plants (Figure 3C).

To get an overview of the expression profile of OsPsbS1, the CREP database (http://crep.ncpgr.cn/crep-cgi/home.pl), a website containing the dynamic gene expression atlas of indica rice, was searched (Wang et al., 2010). It was found that OsPsbS1 exhibited high-level expression in leaf but low-level expression in root and stem in Minghui 63 and Zhenshan 97, two indica cultivars (Figure 4A). In addition, we also investigated the relative expression level of OsPsbS1 in 02428, a japonica variety, by quantitative RT-PCR. Similarly, it was highly expressed in green tissues including leaf and leaf sheath (Figure 4B), indicating OsPsbS1 was a green tissue-specific expressed gene. To complement the results at the mRNA level, a 3141-bp promoter region of OsPsbS1 was ligated in frame with the β-glucuronidase (GUS) reporter gene, which was then introduced into Zhonghua 11, a japonica variety. GUS staining of the P_OsPsbS1::GUS transgenic plants showed that OsPsbS1 had specific and high-level expression in mature leaf (Figure 4E), leaf sheath (Figure 4F), and young glume (Figure 4G), which was in accordance with its high expression in the green tissues of 02428. No staining signals was observed in root (Figure 4C), stem (Figure 4D), and mature seed (Figure 4H). We further compared the expression levels of OsPsbS1 in various tissues before heading stage between two varieties with different genotypes in the OsPsbS1 promoter region by quantitative RT-PCR (Figure 4I). Zhonghua 11 contains the 2,674-bp insertion fragment in the promoter region, whereas Zhenshan 97 doesn't. The results showed that the expression level of OsPsbS1 in Zhonghua 11 was much higher than that in Zhenshan 97 in green tissues such as leaves and shealth (Figure 4I). However, OsPsbS1 exhibited very low-level expression in root and stem in both Zhonghua 11 and Zhenshan 97 (Figure 4I). In conclusion, OsPsbS1 had specific and high-level expression in green tissues of rice, and genotypes with the 2674-bp insertion fragment in the promoter region tend to exhibit much higher level expression.

A candidate gene encoding light-induced protein 1-like (LOC_Os01g01340) was located in one significant association region and close to the representative GWAS peak signal (around 129,845 on chromosome 1) for NPQ_13 in indica group (Table 2). Using resequencing data, the haplotypes of LOC_Os01g01340 were built based on SNPs and InDels detected in the promoter region (around 1 kb upstream of transcription initiation site), 5′ UTR, 3′ UTR, and non-synonymous SNPs and InDels in the coding region. We only focused on those haplotypes shared by at least 5 accessions in the analyzed group. Five haplotypes were found in the whole association panel (Figure S1). In total, five SNPs in the promoter region, three non-synonymous SNPs in the coding region, and seven SNPs in 3′ UTR were detected. Hap1, hap2, and hap3 were the three major haplotypes in indica group, whereas hap1 was the only major haplotype in japonica group (Figure S1). Considering the complexity of population structure and genetic background, we performed statistical analysis in the same or similar groups. Comparisons of the mean values of NPQ_13 showed that hap2 had significantly lower value than hap1 (P = 4.479 × 10⁻¹²) and hap3 (P = 1.227 × 10⁻⁶) in indica group (Figure 5). In addition, hap2 also exhibited significantly lower value than hap1 (P = 0.0019) and hap3 (P = 0.0260) in NPQ_12 in indica. We checked the expression of LOC_Os01g01340 in 59 accessions randomly chosen from indica with different haplotypes. Using qRT-PCR analysis, it was found that expression levels of LOC_Os01g01340 in hap1 accessions were much higher than those in hap2 and hap3 accessions (P < 0.01) (Figure 5, Table S5). These results indicate that this gene might be a good candidate for the GWAS locus.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.