Positional Cloning of the Flowering Time QTL qFT12-1 Reveals the Link Between the Clock Related PRR Homolog With Photoperiodic Response in Soybeans

Flowering time and maturity are important agronomic traits for soybean cultivars to adapt to different latitudes and achieve maximal yield. Genetic studies on genes and quantitative trait loci (QTL) that control flowering time and maturity are extensive. In particular, the molecular bases of E1-E4, E6, E9, E10, and J have been deciphered. For a better understanding of regulation of flowering time gene networks, we need to understand if more molecular factors carrying different biological functions are also involved in the regulation of flowering time in soybeans. We developed a population derived from a cross between a landrace Jilincailihua (male) and a Chinese cultivar Chongnong16 (female). Both parents carry the same genotypes of E1e2E3HaE4 at E1, E2, E3, and E4 loci. Nighty-six individuals of the F2 population were genotyped with Illumina SoySNP8k iSelect BeadChip. A total of 2,407 polymorphic single nucleotide polymorphism (SNP) markers were used to construct a genetic linkage map. One major QTL, qFT12-1, was mapped to an approximately 567-kB region on chromosome 12. Genotyping and phenotyping of recombinant plant whose recombination events were occurring within the QTL region allowed us to narrow down the QTL region to 56.4 kB, in which four genes were annotated. Allelism and association analysis indicated Glyma.12G073900, a PRR7 homolog, is the strongest candidate gene for qFT12-1. The findings of this study disclosed the possible involvement of circadian clock gene in flowering time regulation of soybeans.


INTRODUCTION
Soybeans are an important crop globally and provide quality protein and oil for human and animal consumption (Graham and Vance, 2003). Flowering time and maturity are essential traits in crop cultivation. Soybeans, a typical short-day plant, are highly sensitive to day length. Molecular bases of flowering and maturity genes are important for understanding the mechanism regulating photoperiodic response and adaptation in soybeans (Zhao et al., 2016).
In Arabidopsis, the five major flowering time pathways, e.g., the photoperiod, vernalization, gibberellin, autonomy, and aging, are well characterized. In the photoperiodic pathway, the photoperiodic flowering model of GIGANTEA (GI)-CONSTANS (CO)-FLOWERING LOCUS T (FT) is widely conserved among different plant species. The GI protein is a nuclear protein that links the circadian clock and the flowering time pathway. GI interacting with the Flavin-Binding, Kelch Repeat, F-Box 1 (FKF1) protein can induce CO and FT expression (Shim et al., 2017). In addition, GI can directly or indirectly regulate FT expression by interacting with FT repressors or binding to FT promoter regions in a way independent of CO (Sawa and Kay, 2011). PSEUDO RESPONSE REGULATORS 3,5,7,and 9 (PRR3,PRR5,PRR7,and PRR9), the additional members of the TOC1 family, are core genes participating in additional loops essential for clock function (Mizuno and Nakamichi, 2005). In the soybean genome, PRR family seems to have undergone a different gene expansion process during evolution (Quecini et al., 2007).
The major flowering time gene, E1, encodes a legume specific transcription factor containing a B3-like domain that represses the expression of two florigen genes, GmFT2a and GmFT5a. Four main E1 alleles, e.g., E1, e1-as, e1-fs, and e1-nl, have been identified in soybean cultivars. E1 is the strongest genotype in suppressing the expression of GmFT2a and GmFT5a. The e1-as is a leaky allele that can suppress the expression of GmFT2a and GmFT5a, whereas both e1-fs and e1-nl are non-functional (Xia et al., 2012a). E3 and E4 encode the photoreceptors PHYA3 (Watanabe et al., 2009) and PHYA2 (Liu et al., 2008), respectively, which are involved in the regulation of E1 (Xia et al., 2012a). Intriguingly, the soybean E2 gene was identified to be a GIGANTEA homolog (GmGIa), there are three major haplotypes of E2 (Watanabe et al., 2011). Various allelic combinations at the E1, E2, E3, and E4 loci controlling the absence of or reduced photoperiod sensitivity determine the diversification of soybean maturity and adaptation at different latitudes (Xu et al., 2013;Tsubokura et al., 2013a;Tsubokura et al., 2013b;Zhai et al., 2014;Jiang et al., 2014;Kurasch et al., 2017;Langewisch et al., 2017;Li et al., 2017;Gupta et al., 2017). These genes or loci also influence the yield, adaptability, and quality of soybeans (Xia et al., 2012b;Sayama et al., 2010;Guo et al., 2015). Diversified functions of CO homologs in soybeans have been disclosed; however, these finding lacks the support of QTL study. Intriguingly, soybean E2 gene was identified to be a GIGANTEA homolog (GmGIa). Even now, we still cannot answer to what extent the GI-CO-FT module is conserved in soybean.
In 2018, a QTL for flowering time, qFT12.1, was mapped to a 2,703-kB region interval between the markers BARCSOYSSR_12_0220 and BARCSOYSSR_12_0368 on chromosome 12 (Liu et al., 2018). Although the studies indicated functional genes exist in this region that controls flowering time, it has not been deciphered at high resolution. In this study, a major QTL for flowering time was detected on chromosome 12 in a mapping population derived from a cross between a Chinese landrace and an elite cultivar. Subsequently, fine mapping enabled us to delimit this region to 56.4 kB. Based on allelism and association analysis, Glyma.12G073900, a homolog of the clock gene PRR7, is the strongest candidate gene for qFT12-1, the flowering time QTL. This finding showed the possible involvement of circadian clock gene in flowering time regulation in soybeans.

Plant Material and Phenotyping
An F 2 population of 322 individuals derived from a cross between Jilincailihua (JLCLH) and Changnong 16(CN16) was used for QTL mapping. Both parents are carrying the same genotypes of E1e2E3HaE4 at E1, E2, E3, and E4 loci. Female parent CN16 (E1e2E3HaE4), an elite cultivar that has been widely cultivated in the Northern China, was bred in 2003 at the Changchun Academy of Agricultural Sciences, China. The male parent, JLCLH (E1e2E3HaE4), was a Chinese landrace with a relatively later flowering phenotype. Phenotypic evaluation was conducted at two geographic locations: HARBIN, the campus of Northeast Institute of Geography and Agroecology, Harbin, Heilongjiang (45°70' N, 126° 64' E) and HAILUN, Hailun experiment station of Northeast Institute of Geography and Agroecology, Hailun, Heilongjiang (47° 47N, 126° 97'E). Flowering time (R1) of each plant was defined as days from emergence (DAE) to the opening of the first flower (Fehr et al., 1971). Plots were 5-m long with 0.6-m row spaces. Plants were generally 20 cm apart. We followed standard local agricultural practices to control insects and weeds. Plants were individually harvested.
In 2015, the F 2 population was grown at HARBIN with 178 plants and at HAILUN with 135 plants. Subsequently, subpopulations were phenotyped from 2015 to 2018 at both HARBIN and HAILUN. In 2016, we planted 44 F 2:3 lines whose parent plants had recombination at the QTL region between markers to confirm the authenticity of the QTL detected and to evaluate the suitability for fine mapping the QTL. Based on the analysis of field performance, we selected 4 F 2:3 lines to develop subpopulations for the first round of fine mapping. We further identified seven recombinant plants whose progenies were phenotyped at HARBIN and HAILUN in 2017. Genotyping was performed using 18 markers within qFT12-1 region (Supplementary Table S1). In 2018, the progenies of important recombinants harvested in 2017 were phenotyped for R1 to evaluate the stability of the phenotypic performance.

DNA Isolation
DNA extraction from young leaves was carried out using the CTAB method with minor modifications (Xia et al., 2007). PCR was conducted in a total reaction volume of 10 μl containing 100-ng genomic DNA,10×PCR buffer (with 15 mM/L Mg 2+ ), 0.2 μmol L -1 of each primer, 0.2 mmol L -1 dNTPs, 1 U Taq polymerase (Transgen Biotech, Beijing, China) in distilled water. PCR was performed using the following program: denaturation at 94°C for 5 min, 35 cycles of denaturation at 94°C for 30 s, annealing at 52°C for 30 s, extension at 72°C for 20 s, and a final extension at 72°C for 10 min before cooling down to 16°C.

Genetic Linkage Map Construction and QTL Mapping
The F 2 population was used to construct linkage maps. A total of 96 randomly selected individuals (grown at HAILUN in 2015), along with parents, were genotyped using Illumina SoySNP8k iSelect BeadChip . A total of 7,189 SNPs were specifically manufactured by Infinium HD Ultra. SNP genotyping was performed with the Illumina Iscan platform (Illumina, Inc., San Diego, CA). Genetic linkage analysis and map construction were performed with IciMapping_4.0 (Meng et al., 2015).
The recombination threshold value was set at 0.40 and the Kosambi mapping function was used to convert recombination frequencies into map distances. Markers showing distorted segregation were included in the linkage analysis.
QTL analysis was carried out using the software IciMapping version 4.0. The Inclusive Composite Interval Mapping (ICIM-ADD) model was used to detect QTLs. For the ICIM-ADD, a LOD score of 3.0 was used as a minimum to declare the significance of a QTL in a particular genomic region. One thousand permutations at a probability of 0.05 were also conducted to identify the genome-wide LOD (Meng et al., 2015). QTLs with an LOD score exceeding the genome-wide LOD were declared as significant QTLs, whereas QTLs with LOD less than the genome-wide LOD but more than 3.0 were identified as suggestive QTLs. MCIM was used to map QTLs with additive and epistatic effects. This analysis was performed using a 2D genome scan, with a 1-cM walking speed and 10-cM window size. Significant thresholds (critical F-values) for QTL detection were calculated with 1,000 permutations and a genome-wide error rate of 0.05.
Genotypic and phenotypic data were analyzed using the IBM SPSS Statistics 20 Program. The analysis of variance (ANOVA) of correlation was evaluated under the general linear model in SPSS. Flowering times were tested for deviations from normality using the parameters of kurtosis and skewness by SPSS SPSS 16.0 software (SPSS Inc., Chicago, IL, USA).

Whole-Genome Resequencing of Parents and Marker Development
The parents were sequenced on Illumina HiSeq2000 platform (Illumina, Inc., San Diego, CA, United State), and 125-bp paired-end reads with insert sizes of roughly 300 bp were generated. Forward and reverse paired reads were split for each sample into two separate files (option -split-files). Reads were quality trimmed by NGS QC Toolkit (v2.3.3) (Patel and Jain, 2012) with default parameters. For each sample, only reads passing the quality filtering as matching pairs were retained and aligned to the reference genome Gmax_275_Wm82.a2.v1 using the Burrows-Wheeler Alignment Tool (BWA) v.0.6.2 (Liu and Abe, 2009). The resulting SAM files were converted to sorted BAM files compliant to the Genome Analysis Toolkit (GATK) format by Picard Tools v.1.77 (http://picard.sourceforge. net/) using the tools, in the following order, CleanSam, SamFormatConverter, and AddOrReplaceReadGroups. GATK-compliant BAM files were submitted to GATK v.2.3-3 for preprocessing procedures. The InDel variants above were used to develop Indel markers and SNP variants were used to develop CAPS/dCAPS or High-Resolution Melting (HRM) markers (Supplementary Table S1).
PCR products of InDel primers above were separated using 12% nondenaturing polyacrylamide gels with a 29:1 ratio of acrylamide: bisacrylamide followed by Gelred staining (Transgen Biotech, Beijing, China) and visualized by the GelDoc XR Molecular Imager System (Bio-Rad, USA). HRM primer pairs were designed according to resequence data of parents. HRM ran in the Eco™ Real-Time PCR System (Illumina), and the results were analyzed with Eco Software v3.0.16.0.

RNA Isolation and Quantitative Real-Time PCR (RT-qPCR) Analysis
Total RNA of leaves was extracted as mentioned previously. Reverse transcription was performed with 1 μg of RNA through TransScript ® One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen, Beijing, China) according to the manufacturer's instructions. For the qPCR, 100 ng of total RNA were used for reaction in the Eco™ Real-Time PCR System (Illumina), and the expression levels were analyzed with Eco Software v3.0.16.0 and normalized with the results of soybean TUA5. PCR cycling conditions were set up with the following program: 94°C for 30 s, followed by 40 cycles of 94°C for 5 s, and 60°C for 30 s with fluorescence signal collection. Melting and cooling steps were set as default parameters. Relative gene expression levels were calculated using the 2 −ΔΔCt method. Three independent biological replicates were obtained and subjected to real-time PCR in triplicate. Raw data were standardized as described previously. All primers for expression analysis are listed in Supplementary Table S1.
For the diurnal expression pattern of Glyma.12G073900 and Glyma.U034500, JLCLH and CN16 were grown at LD (16 light/8 dark) and then transferred to continuous light (LL) and dark conditions (DD), respectively. Leaves for RNA extraction were sampled at 3-h intervals at LD and 4-h intervals at DD and LL.

Phylogenetic Analysis
Known PRRs family proteins were collected from the phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) and NCBI (https://www.ncbi.nlm.nih.gov/). The phylogenetic tree was constructed using the neighbor-joining method (Saitou and Nei, 1987). The bootstrap consensus tree inferred from 1,000 replicates is taken to represent the evolutionary history of the taxa analyzed (Felsenstein, 1985). Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. In total, 46 nucleotide sequences including Glyma.12G073900 and its homologues from soybeans and other plants were subjected to phylogenetic analysis. Evolutionary analyses were conducted in MEGA7 (Kumar et al., 2016).

RESULT
Phenotypic Performance for the Parents and the F 2 Population JLCLH, a landrace of semiwild soybean, has the same allelic genotypes of E1e2E3HaE4 as CN16, an improved cultivar, at the E1, E2, E3, and E4 loci (Zhai et al., 2014). However, JLCLH flowered on 65 (days after emergence, DAE), approximately 24 days later than CN16 (flowered on 41 DAE) at HARBIN location (Figures 1A-C). One-way ANOVA indicated that there were significantly difference (p < 0.01) in flowering time and maturity between two parents in both locations (Figures 1D, E). The flowering time in the F 2 population exhibited a broad distribution with a range from 38 to 83 DAE ( Figure 1F). It was hypothesized that there must be other locus or loci influencing flowering time besides E1, E2, E3, and E4.
A Single Major QTL qFT12-1 for Flowering Time Was Mapped to Chromosome 12 Nighty-six individuals of the F 2 population were picked randomly and genotyped with Soy6K_MSU_SNP including 7,189 SNP markers (Akond et al., 2013;Supplementary Table S4). After exclusion of rare allele (a minor allele frequency of less than 10%) or unmapped markers, a total of 2,407 SNP markers were genetically mapped to 20 linkage groups corresponding to the 20 chromosomes of soybeans with QTL IciMapping 4.0 software (Meng et al., 2015) (Figures 2A, B  Basic descriptive statistics of the linkage maps is presented in Table 1. The map spanned 3,425.69 cm, with an average marker interval distance of 1.4 cm. The genetic length of each linkage groups ranged from 79.87 (Chr 3) to 293.23 cm (Chr 12). On average, each linkage group consists of 65 SNP markers, covering 171.28 cm. The recombination rates in the putative centromeric or pericentromeic regions were obviously suppressed in most chromosomes, as shown in the middle portions of each linkage group. In some other regions, gaps also exist due to lacking polymorphic markers between parents, indicating these regions are genetically homogenous between parents.
One QTL was identified for flowering time in the genetic linkage map constructed in the F 2 mapping population with SNP markers. The major QTL qFT12-1 was mapped on chromosome 12 with an LOD score of 4.5, PVE (phenotypic variance explained) of 20.5%, and DOM (dominance effect) of 2.7 ( Table 2). In order to verify the authenticity of qFT12-1, we developed 10 polymorphic markers in the region of qFT12-1. When a total of 171 F 2 plants were used, qFT12-1 was confirmed ( Figure 2C). Moreover, one of markers at nucleotide 5520945 showed significantly correlation with flowering time (P = 0.004), when we analyzed association of genotype and phenotype using 171 F 2 population. We analyzed the two F 2 populations grown in HARBIN (Population A) and HAILUN (Population B) respectively to test the effect of environmental factor on flowering time. Statistical analysis revealed that the environmental factor as well as interaction of environment × genetic significantly affect flowering time (Table 3). Although with a relative lower LOD, qFT12-1 was detected in all environments, revealing qFT12-1 is stable and true QTL. The relative lower LOD might reflect the complexity of flowering time gene network in which many other genes may also play minor or conditional roles in regulation of flowering time. The linkage disequilibrium (LD) was analyzed. The result showed that the LD block is rather large because of the LD (Figure 2D). It is difficult to narrow down the qFT12-1, so we used recombinants to further narrow down qFT12-1.
qFT12-1 Was Delimited to a 56.4-Kb Region From 5496042 to 5552485 on Chromosome 12 To narrow down the genomic position of qFT12-1, the parents were resequenced at a higher coverage level to enhance the chances of developing more SNP and InDel markers. The average sequencing depth of the parents JLCLH and CN16 were 8.1× and 6.5×, respectively. The reference genome Gmax_275_Wm82.a2.v1 comprised 955,380,172 bp of assembled and anchored sequences. After data filtering, 2204598 in JLCLH and 1662895 in CN16 polymorphic SNPs and InDels were high quality, which could be used to design primers. Based on the resequencing data, 64 InDel primers were developed around and between the original map anchor markers Gm12_4670638_A_C and Gm12_5230528_G_A. PCR products amplified from 13 primers pairs showed clear polymorphisms between parents on agarose gel.
In general, the genotyping of 11 recombinants in the region of qFT12-1 was performed for all 18 markers, and the results were indicated in Figure 3. The genotypes of qF12-1 were reflected from phenotypes of progenies of these recombinants, i.e., early flowering homozygous CN16 (qFT12-1), later flowering homozygous JLCLH, and heterozygous JLCLH/ CN16 genotypes. Early flowering phenotypes for the progeny of R121, R121, R103-18, in combination of heterozygous phenotype for R128-27 revealed that the qFT12-1 exist on right side of 24HRM. The late flowering time for progeny of R113-3, early flowering times for progeny of R24-8 and R128-11, and in combination of heterozygous phenotype for R113-23 can delimit qFT12-1 to the left side of 58HRM. Furthermore, the R158 and R171 could narrow down the qFT12-1 region further. From Figure 3, the genotypes of qFT12-1 of R158 and R171 are homozygous JLCLH and homozygous CN16 type judging from the phenotype of progenies of two recombinants, we can delimit the qFT12-1 to a 56.4-kB region between markers of 48HRM and 53 HRM. Furthermore, heterozygous recombinants R113-23 and R128-27 were used to test association between genotype of markers and phenotype in progeny. We genotyped progenies of R113-23(n = 31) and R128-27(n = 32) using 52HRM marker. Homozygous JLCLH and homozygous CN16 genotypes of 52HRM showed significant correlation in flowering time in progenies of recombinants R113-23 and R128-27 (P < 0.03 and P < 0.01, respectively) ( Table 4). Taken together, qFT12-1 was narrowed down to the 56.4-kB region between 48HRM and 53HRM. According to the Williams 82.a2.v1 genome annotation, four genes were annotated: three annotated genes (Glyma.12G073700, Glyma.12G073900, and Glyma.12G074000) having putative functions and Glyma.12G073800 with unknown function ( Table 5).

Glyma.12G073900 Is the Strongest Candidate Gene
When four candidate genes were compared between parents using the resequenced data, there are three SNPs in introns of Glyma.12G073700, one SNP in intron of Glyma.12G073800 and no polymorphic SNP in Glyma.12G074000. However, a total of six SNPs, four in intron and two in exon, have been detected in Glyma.12G073900 (Table 6). Glyma.12G073900 is predicted to encode pseudo-response regulator (PRR) protein with a pseudo receiver domain at the N terminus and a CCT motif at the C terminus ( Figure 4B). PRR genes are clock-associated protein factors known to be key regulator in     photoperiod network in many grains (Murakami et al., 2005;Nakamichi et al., 2010). Among two SPNs occurred in exons of Glyma.12G073900, the first one occurred at the first exon, at nucleotide 490, leading to A in JLCLH to C in CN16. This nonsynonymous mutation caused a change glutamine (Gln) in JLCLH to leucine (Leu) in CN16 at amino acid level. The second occurred at nucleotide 2088, this mutation leading to C in JLCLH to T in CN16, resulting in a premature termination codon (PTC) occurred in CN16 (Figure 4A). We predicted gene interaction on Soynet (https://www.inetbio.org/soynet/), a functional gene network and coexpression network database for a legume crop, soybean (Glycine max). There are mutual regulatory relationships between Glyma.12G073900 and E2, E3, E4, and GmFTs. Glyma.12G073900 and its homologous Glyma.U034500 are interrelated with E2 directly (Supplementary Figure S2). Three-dimensional (3D) structures of supposed PRR proteins from JLCLH and CN16 were predicted. As expected, the putative 3D protein structures of coding protein between parents were different at the C terminus. Putative PRR protein containing a CCT NSC, Non Synonymous coding. motif in JLCLH has four additional beta sheets at C terminus ( Figure 4C) when compared with that in CN16. Therefore, we asserted that Glyma.12G073900 is a candidate gene for qFT12-1 regulating flowering. Phylogenetic tree of PRRs identified in Arabidopsis, rice, indicates that PRR genes Glyma.12G073900 and Glyma. U034500 belong to the same clade ( Figure 5) as PRR3, PRR7 from Arabidopsis. When we compared genetic variations of Glyma.12G073900 in 642 soybean accessions (Supplementary Table S3), 75% of cultivars or accessions are of CN16 genotype (at nucleotide 2088), whereas 23% are of JLCLH genotype (Figure 6). Among 151 cultivar or accession carrying JLCLH genotype, 76 accessions are from Glycine soja, 55 accessions are landraces, and 20 accessions are modern cultivar ( Table 7). In contrast, among 480 CN16-type accessions, 333 accessions are improved cultivars, 126 accessions are landraces, and 21 accessions are from Glycine soja ( Table 7). The result indicated JLCLH allelic type widely exists in Glycine soja while CN16 allelic type is wide distributed in improved cultivars, especially from Northern part of China. 182 of 333 accessions have maturation data, we found 149 accessions are from maturity group 000 to IV, in contrast, only 33 accessions are from maturity group V to X.

Diurnal Circadian Expression Pattern of Glyma.12G073900 and Glyma.U034500
The diurnal expression pattern of Glyma.12G073900 and its homologue in soybeans may help us to clarify whether the candidate gene regulates flowering time through the circadian clock related pathway. Under LD (16 h light/8 h dark), Glyma.12G073900 showed diurnal expression in both JLCLH and CN16, and in both LL and DD condition (Figures 7A-D), the expression peaks were around ZT3 to ZT6 in the parents. Glyma. U034500, the most homologous gene to Glyma.12G073900, showed a similar diurnal expression rhythm, but its expressional abundance is lower than Glyma.12G073900 under LD, and significantly lower under LL and DD.

Genes or QTLs Involved in Regulation of Flowering Time Trait in Soybeans
Flowering time, a vital agronomic trait, influences the productivity, adaptability, and quality of soybeans. Understanding the genetic mechanism of flowering time and photoperiod insensitivity is advantageous for promoting final yield in specific regions. To  Supplementary Table S2. October 2019 | Volume 10 | Article 1303 Frontiers in Plant Science | www.frontiersin.org date, six major E genes for maturity of soybeans (E1, E2, E3, E4, E9, and E10) have been identified. In order to locate novel loci or genes of flowering time in soybeans, parents of the mapping population used in this study were chosen because they had large differences in flowering time although having the same allelic genotypes at E1 to E4. We detected and characterized a major flowering time QTL, qFT12-1, on chromosome 12 in soybeans. Of the 104 flowering time QTLs archived in SoyBase (https://soybase.org/sbt/), one QTL has been reported on chromosome 12 with additive effect of −1.8 to −3.2 days, 26.6-36.8% of PVE at the 19.32-cm interval between Satt568 and Satt442 (Kuroda et al., 2013). Removing the effect of E1, a QTL qFT-H on chromosome 12 with LOD 2.38 with nearest marker of Satt442 in 28.2-cm interval was detected (Lu et al., 2016). Kuroda et al. (2013) identified the flowering time QTL including the region of qFT12-1 by using F 2 populations derived from crosses between a wild soybean accession (JP110755) and an improved soybean cultivar, Fukuyutaka. Liu et al. (2018) mapped qFT12.1 to a 2703-kB (44.04 cm) interval between the markers BARCSOYSSR_12_0220 and BARCSOYSSR_12_0368 on chromosome 12 with an LOD score of 15.16 and PVE of 38.27% using a CSSL population derived from a soybean cultivar, Jackson (PI548657) and a wild soybean accession, JWS156-1. The existence of QTL for flowering time frequently detected in different populations and environments suggests that a key gene in this region on chromosome 12 plays important role in controlling flowering time. Besides flowering time QTLs, two more pod maturity QTLs were detected around this region. Moreover, a genome-wide association study (GWAS) showed that seven more flowering time QTLs on chromosome 12 .

qFT12-1 Region Has Been Subjected to Extensive Selection During Evolution and Domestication
Modern cultivated soybean was widely believed to domesticate from wild soybean (Glycine soja) (Kim et al., 2012). It has been proven that wild soybean or landraces preserve higher genetic diversity including genes controlling flowering time, providing genetic resource pools for the improvement of cultivated soybeans   (Zhou et al., 2015). Domestication could increase plant adaptability to different environments (Allaby, 2010). In sorghum, there are many loss-of-function alleles encoding truncated PRR37 protein.
Each PRR37 allelic variant could be traced to specific geographic location or specialized agronomic type (Klein et al., 2015), indicative of human selection. In this study, when we analyzed the variation of Glyma.12G073900 at nucleotide 2088 in 642 soybean accessions, it was found that JLCLH homologous type mostly exists in wild soybeans and landrace, whereas CN homozygous type are mostly distributed in improved cultivar. An inDel located on 5520259 (reference genome position) in Glyma.12G073900 has been identified, which can clearly distinguish wild soybeans from modern soybean cultivars when whole genome resequencing of another batch of soybean germplasms were analyzed (Song QJ, personal communication). This result reassure that this gene is related to domestication and the region is within domestication sweep. Among 182 soybean accessions or cultivars being assigned to the maturity groups (Supplementary Table S3), 149 accessions or cultivars were belonged to MG 000 to MG IV suitable for cultivation in high latitude. The allelic type Glyma.12G073900 in JLCLH is probably ancestral genotype and this gene has been subjected to natural selection in evolution and artificial selection during domestication process. According to the analysis of resequencing data of 302 wild and cultivated accession of soybean, the region of chromosome12: 5271201-58199864 is within domestication sweep, in the proximal region of qFT12-1, a gene controlling pubescence has been subjected to intensively artificial selection. The presence of qFT12-1 have been reported by different researchers using different populations all of which were derived wild soybean × improved cultivars (Kuroda et al., 2013;Liu et al., 2018). We assumed that CN16 homozygous genotype of qFT12-1is generally required for modern cultivars to reduce photoperiod sensitivity to adapt to higher latitudes.

PPR Homolog May Be Involved in Circadian Clock and Flowering Time Regulation in Soybean
PRRs are key components of transcription/translation circadian networks in plants (Farré and Liu, 2013). PRRs are core genes in the translational feedback loop mechanism of the central oscillator in Arabidopsis (Hayama and Coupland, 2003;Nakamichi et al., 2007). PRR9, PRR7, and PRR5 of Arabidopsis were identified to suppress CCA1 and LHY promoter activities (Nakamichi et al., 2010). PRR proteins with a pseudo receiver domain at the N terminus and a CCT motif at the C terminus were generally proposed to contribute to day-length measurement through control of the time-keeping mechanism associated with CO transcription (Yamamoto et al., 2003;Murakami et al., 2004). OsPRR37 was mapped very closely to Hd2-QTL, which was identified as the major locus that enhances the photoperiod sensitivity of flowering (Murakami et al., 2005). In sorghum accession ATx623, one additional nucleotide substitution in SbPRR37 resulted in a truncated protein without CCT motif, which leaded to photoperiod insensitive and early flowering (Murphy et al., 2011). Expression of floral repressors SbPRR37 flowering under LD conditions inhibits flowering in sorghum . In this study, the candidate gene for qFT12-1, Glyma.12G073900 is the homogenous to PRR7 of Arabidopsis. A single mutation occurred at the end of the eighth exon, Glyma.12G073900 in CN16 resulting in a CCT motif truncated protein. CCT motif might be the key domain for its function in regulation of flowering. Circadian rhythms in expression of Glyma.12G073900 both in JLCLH and CN16 are well accordance with that PRR domains in Glyma.12G073900 are conserved in JLCLH and CN16. As short day plant, soybean has a special gene network for regulation of flowering which is different from Arabidopsis. In soybean, E1 is the most important gene in flowering regulation network. In flowering pathway of soybean, E1 accepts light signal from phytochrome through some circadian clock genes, regulates flowering by expression of GmFTs (GmFT2a, GmFT5a, GmFT4). In this study, the candidate gene of qFT12-1, Glyma.12G073900 has characteristic of circadian clock gene and interactively interoperate with GmFTs as showed prediction network. If it had, Glyma.12G073900 most probably conducts its function with E2 encoding GIGANTEA. More research is needed to verify this hypothesis and how qFT12-1 coordinately controls flowering time in cooperation with the E1 pathway.

CONCLUSION
In the present study, a major QTL for flowering time, qFT12-1, was identified, validated, and characterized. At first, we mapped the qFT12-1 to an approximately 567-kB region on chromosome 12. Using recombinants whose recombination events were occurring within the QTL region allowed us to narrow down the QTL region to 56.4 kB. The PRR7 homolog Glyma.12G073900 is proven to be the strongest candidate gene for qFT12-1. Further studies are needed to determine the effect of qFT12-1 in other genetic backgrounds, as well as its interaction with other flowering time genes and with environmental conditions, such as temperature and light quality. Regardless, the deciphering of molecular basis of qFT12-1, a QTL of flowering time will not only contribute to a better understanding of the molecular basis of the photoperiod-dependent regulation of flowering time but could also be potentially useful in improving the adaptability and productivity of cultivated soybeans.

ACKNOWLEDGMENTS
Thanks to Dr Zixiang Wen from Michigan State University for the help in genotyping using Illumina SoySNP8k iSelect BeadChip.