Minor Effects of 11 Dof Family Genes Contribute to the Missing Heritability of Heading Date in Rice (Oryza sativa L.)

DNA binding with one finger (Dof) proteins are plant-specific transcription factors with important and diverse functions in seed germination, flowering time, and biotic and abiotic stresses. In this study, haplotype-based association analysis was conducted between heading date and 30 Dof family genes in a worldwide germplasm collection. Of these, 22 Dof genes were associated with heading date. Multiple comparisons among haplotypes revealed their diverse functions in promoting and suppressing heading date under short-day (SD) and long-day (LD) conditions. They cumulatively made a considerable contribution to the missing heritability of heading date. A set of knockout mutants of 30 Dof genes generated by CRISPR/Cas9-mediated genome editing technology showed that 11 and 9 Dof genes regulated heading date under LD and SD, respectively. Phenotype measurement of mutants showed that these 11 and 9 Dof genes slightly regulated heading with effects of 2–5 days under LD and SD, respectively. Both mutant and natural variation assays indicated functional redundancy in regulating heading date among Dof family genes. Nucleotide diversity analysis suggested that most Dof genes have been subjected to selection during domestication and improvement. Beyond heading date, this set of mutants is also a good resource for evaluating the function of Dof genes in regulating stress tolerance and seed germination.

Functional characterization of gene families is often made by overexpression, knockout (gene editing), or knockdown (RNA interference). Overexpression sometimes causes ectopic expression and does not necessarily reflect the native function of target genes. Knocking out and knocking down genes are better ways to identify gene function because these methods precisely target the gene in situ, which does not cause genetic noise. However, the effect of manipulating single genes is reduced if functional redundancy exists within the gene family. Currently, the clustered regularly interspaced short palindromic repeats (CRISPR) system is a powerful tool to produce mutants in target genes (Cong et al., 2013;Hwang et al., 2013;Shan et al., 2013;Yin et al., 2017;Zhou et al., 2018). CRISPR technology relies on Cas9 and single guide RNA (sgRNA) to achieve precise cutting. CRISPR/Cas9 systems have been applied to multiplexed genome editing Zhang et al., 2016), which will facilitate functional analysis of gene families. Transgene-free plants can be achieved through backcrossing and used directly in breeding. For example, the non-transgenic, low-gluten wheat lines were obtained through modification of the a-gliadin family genes using CRISPR/Cas9 technology and serve as source material to introgress this trait into elite wheat varieties (Sanchez-Leon et al., 2018).
Next-generation sequencing technology has developed rapidly and has greatly facilitated genome-wide association studies (GWAS) (Huang et al., 2010;Huang et al., 2017), which have been widely used to identify candidate genes associated with various complex traits in several species including humans, animals, and plants (Wellcome Trust Case Control C, 2007;Moreira et al., 2018;Su et al., 2018). Candidate gene-based association analysis is an option to identify potential functions of related genes. Moreover, association mapping at the haplotype level is a very effective method to establish the relationships between genes and traits Chen et al., 2018;Lu et al., 2018). Haplotype development is easy at the single gene level, which is helpful to identify the association between related genes and traits at the haplotype level.
DNA binding with one finger (Dof) proteins are plantspecific transcription factors (Yamagishi et al., 2002;Yamaguchi et al., 2004). Since the first Dof gene, MNB1, was identified in maize (Yanagisawa and Izui, 1993), an increasing number of Dof genes have been identified in various plant species. Recent studies have revealed the role of the Dof family genes in multiple plant developmental processes. In tomato, TDDF1 (TOMATO DOF DAILY FLUCTUATIONS 1) is involved in circadian regulation and stress resistance (Ewas et al., 2017). In maize, ZmDof3 regulates starch accumulation and aleurone development in maize endosperm (Qi et al., 2017). In rice, a total of 30 Dof genes have been predicted through genome analysis (Lijavetzky et al., 2003). All have a highly conserved DNA-binding domain (Dof domain) of 50 amino acid residues including a C2C2-type zinc finger motif (Huang et al., 2018). RDD1/OsDof2 overexpression promotes flowering in addition to regulating nutrient ion uptake and accumulation (Iwamoto et al., 2009;Iwamoto and Tagiri, 2016). OsDof4 has distinct flowering effects under long-day and short-day conditions . OsDof12 participates in the regulation of flowering time and plant architecture (Li et al., 2009;Wu et al., 2015). These studies indicate that Dof family genes in rice might have important functions in flowering. Therefore, we hypothesized that some other Dof family genes are probably related to flowering. To test this hypothesis, association analysis between heading date and haplotypes of Dof family genes were conducted in this study. On the one hand, candidate gene-based association analysis of Dof family genes with heading date was examined at the haplotype level in a worldwide germplasm collection of 529 rice accessions. On the other hand, the CRISPR/Cas9 gene editing technology was utilized to test which Dof family genes really function in the regulation of heading date. Our results showed that 11 and 9 Dof genes regulated heading date under long day (LD) and short day (SD), respectively. Nucleotide diversity analysis showed that most Dof genes regulating heading date have been subjected to selection during domestication and genetic improvement.

Plant Materials and Field Experiments
Zhonghua 11 (ZH11; Oryza sativa japonica) was used as the recipient for transformation. A total of 529 rice accessions from around the world, composed of a core Chinese collection of 202 cultivars and a core world collection of 327 cultivars, were grown in Wuhan (114°21′E, 30°28′N; the average day length is more than 13.5 h from the middle of May to the beginning of August; 23-31°C) in the summer of 2013 and Hainan (110°01′E, 18°30′ N; the average day length was less than 12.5 h from December to the middle of March; 20-28°C) in the winter of 2013 (https:// www.timeanddate.com, https://www.climatestotravel.com). The details of field management and measurement of heading date are described in a previous study . The 529 accessions were genotyped via sequencing and classified into four subpopulations: indica, aus, japonica, and Admixture (Chen et al., 2014). Single nucleotide polymorphism (SNP) information about the 529 cultivars is available on RiceVarMap v1.0 (http:// ricevarmap.ncpgr.cn/v1/), which is a comprehensive database of rice genomic variations. All transgenic plants were grown in the field in Wuhan (LD) in summer 2018 and Hainan (SD) in spring 2019. Heading date was individually scored as the number of days from sowing to the emergence of the first panicle on the plant. Plant height was measured from the surface to the top of the main panicle.

Haplotype Construction and Haplotype-Based Association Analysis
To obtain sequencing data for all Dof genes in rice, the conserved Dof domain sequence of the known protein SP3 was used to search the Rice Genome Annotation Project (http://rice. plantbiology.msu.edu) and NCBI (https://www.ncbi.nlm.nih. gov/) (Huang et al., 2018). A total of 30 Dof genes were identified.
For development of haplotypes, first, SNPs in Dof genes were extracted from the sequences containing a 2-kb promoter region and the gene body region. Then, PHASE software, which can estimate missing and heterozygous genotypes, was used to construct haplotypes (Stephens et al., 2001). Finally, the haplotypes with an allele frequency ≥0.01 (five accessions) were included for association analysis. Analysis of variance (ANOVA) was used to identify the association between haplotype and heading date. Duncan's test was employed to perform multiple comparisons between all possible haplotype pairs. The phenotypic variation explained by these 22 Dof genes was calculated by ANOVA with a linear model including all these genes.

Nucleotide Diversity Analyses
DNA sequences of each Dof family gene in 1,612 cultivar accessions and 446 Oryza rufipogon accessions were obtained from the ECOGEMS database (http://ecogems.ncpgr.cn). Each sequence contained a 2-kb promoter region and a gene body region. The details of these accessions and their sequencing data have been previously reported (Huang et al., 2012a). The nucleotide diversity (p), Tajima's D, and Fu and Li's D statistics were calculated both in cultivars and wild rice (O. rufipogon) using the DnaSP 6.0 program (Rozas et al., 2017).

Vector Construction and Transformation
To generate knockout mutants of Dof family genes by CRISPR/ Cas9, 19-and 20-bp fragments in the exon regions of Dof genes were chosen as the candidate targets according to the design principles of the target sequences in the CRISPR/Cas9 system (Supplementary Table S1). At least one target was in the exon regions of each Dof gene (Supplementary Figure S1). The two target fragments were introduced into two sgRNA expression cassettes by dual-nested PCR, driven by the OsU6a and OsU6b promoters, respectively. Then, the multiple sgRNA expression cassettes were ligated to the CRISPR/Cas9 binary vector (pYLCRISPR/Cas9Pubi-H) based on Golden Gate cloning . Finally, 30 recombinant CRISPR/Cas9 constructs were introduced into Agrobacterium tumefaciens strain EHA105 and separately transferred into ZH11 by Agrobacteriummediated transformation (Lin and Zhang, 2005).

Mutation Detection and Transgene-Free Line Screening
T 0 transgenic plants were used to detect mutations. Genomic DNA was extracted from the leaves of T 0 plants using the cetyltrimethyl ammonium bromide (CTAB) method. Primer pairs flanking the designated target sites were used to amplify the potentially mutated fragments. PCR products were sequenced to detect mutations. Homozygous mutations were identified using the Sequencher 5.1 Demo software and Degenerate Sequence Decoding method (DSDecodeM; Liu et al., 2015). In our study, the homozygous transgenic knockout lines were transgene-free (or transgene-clean) mutant lines, which were obtained by screening the recombinants between Cas9 and the targeted mutant gene in the selfing progeny. Thus, the phenotypic evaluation was performed between transgene-clean mutant lines and wild type, not the negative plants. The seeds were harvested from each homozygous T 0 plant, and the transgene-free plants were identified from T 1 homozygous families by agarose gel electrophoresis to separate PCR products with Cas9-specific primers. Cas9-negative plants were transgene-free. CRISPR/ Cas9 plasmid DNA was used as a positive control, and ZH11 DNA and H 2 O were used as negative controls.

Gene Structure and Haplotype Analysis of Dof Family Genes
The 30 Dof genes were not evenly distributed in the genome. Chromosomes 1 and 3 each had the largest number of genes (six), and no Dof genes were located on chromosome 11 ( Table 1). In general, Dof family genes had relatively simple structures. Of all Dof genes, 17, 12, and 1 had no intron, one intron, and two introns, respectively (Supplementary Figure S1). The number of SNPs in the 2-kb promoter and entire genomic sequence of Dof family genes ranged from 14 (OsDof14) to 113 (OsDof10 ; Table 1). Based on these SNPs, we constructed the haplotypes of all 30 genes in the germplasm collection of 529 accessions. In the full population, the number of haplotypes ranged from 4 to 20, among which OsDof10 had the largest number (twenty) of haplotypes, whereas OsDof18 had the lowest number (four) of haplotypes. Additionally, with the exception of OsDof4, OsDof21, OsDof29, and OsDof30, all Dof genes had more haplotypes in the indica subpopulation than the japonica subpopulation ( Table 1).

Comparison of Haplotype Effects on Heading Date
Multiple comparisons of haplotypes were made for heading dateassociated genes. Under LD, significant differences in heading date were detected among major haplotypes of 10 genes (OsDof1,6,7,10,11,12,13,16,20,and 21) in both indica and japonica accessions (Supplementary Table S2). Under SD, a similar result was also observed in these 10 genes except for OsDof1 and OsDof12 (Supplementary Table S2).
OsDof12 had six major haplotypes in the germplasm collection. Hap1, Hap2, and Hap3 were mainly harbored by the indica subpopulation, Hap4 and Hap5 were primarily harbored by the japonica subpopulation, and most aus accessions belonged to Hap6 (Figure 2A). Under LD, Hap2 (104.4 ± 5.9 days) and Hap3 (101.0 ± 10.7 days) significantly delayed heading date compared to Hap1 (91.1 ± 9.1 days) within indica rice ( Figure 2B), and Hap4 (114.8 ± 10.9 days) significantly delayed heading date compared to Hap5 (94.2 ± 13.2 days) in japonica rice ( Figure 2C). These functionally diverse haplotypes provide us with an opportunity to improve rice heading date for specific ecotypes. However, under SD, there were no significant differences in heading date among major haplotypes in either indica or japonica rice ( Figures  2B, C). Interestingly, only one SNP (sf0303739795) causing a non-synonymous mutation (G267A) was detected in the coding region of OsDof12, and most aus, indica, and japonica accessions carried Gly at site 267 (Hap2-Hap6). Only a small minority of indica accessions carried Hap1 with Ala at site 267 ( Figure 2A). Within indica rice, Hap2 with G267 and Hap3 with G267 significantly delayed heading date compared to Hap1 with A267 in LD ( Figure 2B). Therefore, G267A might be the causal mutation for the variation of heading date in indica rice. Accordingly, OsDof1 had significant effects on heading date in both indica and in japonica rice under LD, whereas no significant effects were observed under SD (Supplementary Table S2).
In addition, five genes (OsDof2, 23, 25, 29, and 30) had significant effects in japonica accessions and not in indica rice under SD (Supplementary Table S2). These five genes, except OsDof2, also had similar effects under LD (Supplementary Table S2). In the case of OsDof23, six major haplotypes were obtained. Hap1, Hap2, and Hap3 were mainly harbored by the indica subpopulation, whereas Hap4, Hap5, and Hap6 were primarily harbored by the japonica subpopulation ( Figure 3A). Two SNPs (sf0729076798 and sf0729080067) in the coding region caused amino acid changes. The SNP at site sf0729076798 caused a non-synonymous mutation (D400N). Most aus, indica, and japonica accessions were divided into Hap1-Hap5 and shared an Asp at site 400. A small proportion of japonica accessions possessed Hap6 with Asn at site 400 ( Figure 3A). No significant effects were detected in indica rice in either LD or SD ( Figure 3B), whereas the effects were significant in japonica rice ( Figure 3C). Hap6 with N400 (69.2 ± 9.9 days under LD and 81.6 ± 4.3 days under SD) significantly promoted heading date compared to Hap4 with D400 (100.4 ± 5.4 days under LD and 105.5 ± 7.2 days under SD) and Hap5 with D400 (104.8 ± 13.5 days under LD and 89.9 ± 7.7 days under SD; Figure 3C). Therefore, sf0729076798 might be a functional nucleotide polymorphism site in japonica rice.

Mutagenesis of Dof Genes by CRISPR/ Cas9
Multiple positive plants for each gene were harvested from T 0 transgenic plants of all Dof family genes generated by the CRISPR/Cas9 system (Supplementary Table S3). Transgenefree homozygous knockout lines with various mutations in 30 Dof genes were identified by sequencing. The mutations included small fragment insertions (1 bp), deletions (1-17 bp), and large fragment deletions (111-2,371 bp; Supplementary Figure S3). Because the two target sites were designed in the exon region (at least one target fell in the coding region) of Dof genes, each mutation occurred in the coding region and led to a frameshift or truncation mutation in the protein that may completely deactivate the protein function.

Identification of Dof Family Genes Regulating Heading Date Using Mutants
At least three independent mutant lines for each Dof gene were analyzed for heading date. Compared with the wild-type ZH11, the CRISPR/Cas9 mutated plants of four Dof family genes (OsDof1, 11, 21, and 29) promoted heading by approximately 3-4 days under LD ( Figures 4A, B, E, Table 3, and Supplementary Figure S4). Under SD, a similar result was also observed in these four genes except for OsDof11 and OsDof29 (Supplementary Table S4 and Figure S4). Mutants of seven genes (OsDof2,8,9,16,22,24,and 26) delayed heading by 2-5 days under both LD and SD (Figures 4C-E, Table 3,  Supplementary Table S4, and Figure S4). The Dof family gene mutants also had effects on plant height (−14.6 to +14.1 cm under LD; −22 to +4.9 cm under SD; Figure 4F, Table 3, and Supplementary Table S4). Other members did not show any difference in heading date compared to the wild type (Table 3).

Nucleotide Diversity of Dof Genes
We analyzed the nucleotide diversity of Dof family genes in 446 wild rice and 1,612 cultivars. There were 2,822 SNPs across the entire sequence (including 2 kb of promoter and the gene body region) of the Dof family genes among 2,058 accessions (Supplementary Table S5). The pairwise nucleotide diversity ranged from 7 to 59 SNPs per kilobase for Dof genes (Supplementary Table S5). The average nucleotide diversity (p = 2.4 × 10 −3 ) of Dof family genes in both O. sativa (p = 1.9 × 10 −3 ) and wild rice (p = 2.9 × 10 −3 ) was equivalent to that of the full O. sativa genome (p = 2.4 × 10 −3 ). The ratio of nucleotide diversity in cultivars (p c ) to nucleotide diversity in wild rice (p w ) for OsDof2, 4, 10, 16, 18, and 24 was less than 0.5    , 8, 12, 14, 17, 18, 22, and 29 in cultivars and OsDof3, 5, and 7 in wild rice deviated significantly from neutrality (P < 0.05) (Supplementary Table S5). These results suggested that these genes have been subjected to natural selection in wild rice and artificial selection in cultivar rice.

Minor Effects of Dof Genes on Heading Date Has Contributed to the Missing Heritability in Rice
According to the "omnigenic" model, a complex quantitative trait is controlled by a few core genes (major genes) and a large number of peripheral genes (minor genes), and the sum of minor effects across peripheral genes can far exceed the genetic contribution of core genes (Boyle et al., 2017). This statement is supported by the so-called missing heritability (Manolio et al., 2009). Much missing heritability is due to peripheral genes whose effects are too small to reach the level of genome-wide significance (Yang et al., 2015). Heading date shows high heritability in linkage mapping populations and germplasm collections. Moreover, these core heading date genes explain only a fraction of the total heritability in both GWAS and biparental mapping populations (Xue et al., 2008;Huang et al., 2010;Yan et al., 2011;Huang et al., 2012b;Yan et al., 2013;Han et al., 2016;Yano et al., 2016). There must be many genes with minor effects spread widely across the genome. It is often difficult to identify unknown genes with minor effects due to weak association signals in GWAS at the SNP level and stringent thresholds in linkage analysis. Like our Dof family genes, whose effects are too small to be significantly associated with heading date in GWAS (Supplementary Figure S2). Heading date shows high heritability (0.75), but quantitative trait loci (QTLs) revealed by GWAS only explain 40.6% of the phenotype variation . In this study, 22 Dof genes were associated with heading date via candidate genebased association analysis at the haplotype level (Supplementary Figure S4). They all together explained about 18.4% of the variation in heading date. Moreover, mutants in 11 Dof genes showed a 2-to 5-day difference in heading date compared to wild-type ZH11. It is very likely that there is gene redundancy among Dof family genes, which results in false negatives. In addition, it is possible that some Dof genes in ZH11 (the transformation background genotype) are weak-functional or non-functional because ZH11 is not strongly sensitive to photoperiod (Xue et al., 2008;Yan et al., 2013;Zhang et al., 2015a). Therefore, it is understandable that some knockout mutants show no phenotype change or minor effect. But at least 11 Dof genes probably regulate heading date, none of which has been previously mapped, with the exception of OsDof2. These Dof genes with minor effects adequately explain part of the missing heritability. In fact, there are CCT family genes and HAP family genes in rice, and some of them also control heading date, but they have not been mapped by QTL analysis (Zhang et al., 2015b;Li et al., 2016). These minor genes together contribute to the missing heritability of heading date.

Dof Family Genes Have Diverse Responses to Long-Day and Short-Day Conditions
In this study, we detected many Dof genes related to heading date using haplotype-based association analysis under LD (Wuhan) and SD (Hainan). Of them, 16 genes (OsDof3,6,7,8,10,11,12,13,16,20,21,22,23,25,29,and 30) in the full population, four genes (OsDof6,11,16,and 20) in indica, and six genes (OsDof7,10,11,20,23,and 25) in japonica were commonly detected under both LD and SD (Table 2), which indicated that these Dof genes simultaneously contributed to the heritability of heading date under both conditions. One gene (OsDof1) in the full population, six genes (OsDof1, 3, 7, 8, 13, and 21) in indica, and one gene (OsDof12) in japonica were detected under LD ( Table 2), suggesting that these genes primarily function under LD. These results were supported by a previous report on OsDof12 in which overexpressed OsDof12 promoted rice flowering under LD, whereas overexpression of OsDof12 had no effect under SD (Li et al., 2009). In addition, five genes (OsDof2,14,19,26,and 27) in the full population, three genes (OsDof2, 10, and 14) in indica, and eight genes (OsDof2,6,19,21,22,27,29,and 30) in japonica were detected under SD (Table 2), which implied that these genes primarily function under SD. Furthermore, OsDof2 lossof-function mutants generated by CRISPR/Cas9 showed a delay in flowering time compared with wild-type plants under both LD and SD, which was consistent with previous studies (Iwamoto et al., 2009;Iwamoto and Tagiri, 2016). However, OsDof2 was only detected under SD in this study ( Table 2). Taken together, Dof family genes have distinct responses to photoperiod.

There Is Functional Redundancy and Differentiation Among Dof Family Genes
In general, functional redundancy has been reported among related genes in rice (Zhang et al., 2015b;Li et al., 2016). When knockdown or knockout mutagenesis is used for studying the function of related genes, gene redundancy can cause unclear results. In our study, there were 19 genes ( Table 2) whose knockout mutants had heading date that were indistinguishable from wild type. Phylogenetic analysis showed that six pairs of genes, OsDof1 and OsDof29, OsDof2 and OsDof23, OsDof9 and OsDof18, OsDof12 and OsDof26, OsDof13 and OsDof30, and OsDof24 and OsDof25, share a very conserved gene structure outside the Dof domain (Huang et al., 2018), which is probably responsible for their similar functions. For example, knockout mutants of both OsDof1 and OsDof29 promoted flowering time under LD (Table 3). Additionally, some Dof genes may have functional redundancy, which causes single mutants to have wild-type phenotypes. For example, knockout mutants of both OsDof13 and OsDof30 did not show any change in heading date from the wild type under LD. To specify which genes are functionally redundant, double or triple mutants are needed for phenotype analysis. In addition, knockout mutants of OsDof2, OsDof9, OsDof24, and OsDof26 all show delayed heading under LD, whereas mutants in their corresponding evolutionarily most closely related genes OsDof23, OsDof18, OsDof25, and OsDof12 did not show any difference in heading date, suggesting a possible functional differentiation between them. It has been reported that some Dof genes are probably involved in stress resistance (Corrales et al., 2014). Therefore, the mutants generated in this study can also be used to evaluate stress resistance.

CONCLUSION
In summary, out of 30 Dof genes, haplotype-based association analysis identified 22 genes associated with heading date either under long-day or short-day conditions. Favorable alleles were provided for improvement of varieties with different ecotypes. A set of knockout mutants of all 30 Dof genes showed that 11 and 9 Dof genes regulated heading date under long-day and short-day conditions, respectively. Of them, nine genes regulated heading date under both conditions. Functional redundancy might mask the phenotype of single mutants and result in missing identification of flowering Dof genes. The minor effects of Dof family genes that could not be mapped using GWAS make considerable contribution to the missing heritability of heading date in rice.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
YX and YH designed the experiments. YH performed most of the experiments. YH and ZH analyzed and summed the data and drew the figures. NC, ML, and XB assisted in experiments and discussed the results. YH wrote the manuscript. YX revised the manuscript.