QTL-seq analysis identifies two genomic regions determining the heading date of foxtail millet, Setaria italica (L.) P.Beauv.

Heading date is an important event to ensure successful seed production. Although foxtail millet (Setaria italica (L.) P.Beauv.) is an important foodstuff in semiarid regions around the world, the genetic basis determining heading date is unclear. To identify genomic regions regulating days to heading (DTH), we conducted a QTL-seq analysis based on combining whole-genome re-sequencing and bulked-segregant analysis of an F2 population derived from crosses between the middle-heading cultivar Shinanotsubuhime and the early-heading cultivar Yuikogane. Under field conditions, transgressive segregation of DTH toward late heading was observed in the F2 population. We made three types of bulk samples: Y-bulk (early-heading), S-bulk (late-heading) and L-bulk (extremely late-heading). By genome-wide comparison of SNPs in the Y-bulk vs. the S-bulk and the Y-bulk vs. the L-bulk, we identified two QTLs associated with DTH. The first QTL, qDTH2, was detected on chromosome 2 from the Y-bulk and S-bulk comparison. The second QTL, qDTH7, was detected on chromosome 7 from the Y-bulk and L-bulk comparison. The Shinanotsubuhime allele for qDTH2 caused late heading in the F2 population, whereas the Yuikogane allele for qDTH7 led to extremely late heading. These results suggest that allelic differences in both qDTH2 and qDTH7 determine regional adaptability in S. italica.

Introduction

Although millets are inconspicuous in comparison with major cereals including rice, maize and wheat, millets provide important food source worldwide. They are adapted to a wide range of conditions such as those found in arid, hot, salty or low nutrient environments (Goron and Raizada 2015). Foxtail millet (Setaria italica (L.) P.Beauv.), a diploid grass species (2n = 18) with a relatively small genome size (~515 Mb), is one of the oldest cereals in the world. Because foxtail millet is rich in genetic diversity (~9000 varieties), this species has been widely adapted to various environmental regions particularly in Asia and Africa (Reddy ) and has become an important foodstuff in semiarid regions as found in China and India (Zohary ). Recently, reference genome sequences were published (Bennetzen , Zhang ) that allow us to use genetic approaches for improvement of its agronomic traits.

Heading time is an important developmental transition in plants leading to successful sexual reproduction and is determined by multiple genes and environmental factors, such as day-length and temperature. Plants are classified into three types, long-day plants, short-day plants and day-neutral plants according to their photoperiodic flowering responses. Foxtail millet is short-day plants (Thomas and Vince-Prue 1996). Takei and Sakamoto (1987, 1989) investigated heading traits of foxtail millet landraces collected from various regions of Europe and Asia, ranging from low latitudes, such as Indonesia and other Southeast Asian Countries, to high latitudes, such as Belgium and Kirghizia, and found large variability of heading date in foxtail millet. They also classified the landraces into three types (Type I–III) based on combinations of length of basic vegetative growth and sensitivity to day-length. Although heading is an important factor in the adaptation of crops to a cultivation area, the molecular mechanisms underlining foxtail millet heading remain unknown. Understanding heading mechanisms of foxtail millet is important since it would allow us to modify heading date that may potentially increase both grain yield and improve grain quality.

Forward genetic approaches such as linkage and association mapping have been used in the genetic analysis of foxtail millet. Although many quantitative trait loci (QTL) have been identified for agronomic traits such as heading date, biomass, spikelet-tipped bristles and yield (Fang , Jia , Mauro-Herrera , Mauro-Herrera and Doust 2016, Sato , Zhang ), the QTL intervals have been often large. Jia identified three QTLs for the trait of days to heading (DTH) in foxtail millet using genome-wide association studies and found genes orthologous with rice Heading date 1 (Hd1) (Ishikawa , Yano ), Oryza sativa Pseudo-Response Regulator 37 (OsPRR37) (Koo ) and FLAVIN-BINDING, KELCH REPEAT, F-BOX 1 (FKF1) (Han , Matsubara ) located in each candidate region. Zhang have been identified two QTLs for DTH in foxtail millet and one of two QTLs corresponding to Hd1. Hd1 gene of foxtail millet that is related to latitudinal adaptation and domestication (Fukunaga , Liu ). Investigation of a recombinant inbred line (RIL) population derived from a cross between foxtail millet and green millet (Setaria viridis) identified 16 QTLs for flowering time (Mauro-Herrera ).

Although conventional QTL mapping has been the primary way to detect QTLs for interesting traits, this method requires genotyping of many individuals in a segregating population using a large number of DNA markers. Bulked segment analysis (BSA) was originally developed to rapidly map genes in segregating populations (Michelmore ). QTL-seq, a new technique that combines the BSA method with next-generation sequencing (NGS), can rapidly identify genomic regions containing QTL without a requirement for DNA marker development (Takagi ). The QTL-seq method has been used to detect QTLs in rice, cucumber, tomato, groundnut and foxtail millet (Illa-Berenguer , Lu , Masumoto , Pandey , Wei ).

In this study, we used the QTL-seq method to identify QTLs that regulate heading in foxtail millet. To understand genetic basis for heading in foxtail millet, we used a segregating population derived from a cross between two Japanese cultivars, Yuikogane and Shinanotsubuhime. These cultivars were chosen since the cropping area for Yuikogane is more northward than that for Shinanotsubuhime, indicating that their heading dates determine regional adaptability of the two cultivars. Variation between these two cultivars (Yuikogane and Shinanotsubuhime) is not so large in the entire variation of heading date among a word-wide foxtail millet landrace collection. We report here the chromosomal locations and genetic effects of the QTLs associated with the difference in DTH between Yuikogane and Shinanotsubuhime.

Materials and Methods

Plant materials

Two foxtail millet cultivar lines, Yuikogane and Shinanotsubuhime, were the parental lines used to develop an F2 population segregating for heading date. Yuikogane was bred and released as a new foxtail millet cultivar in Iwate prefecture, a northern region of Japan (39°42′13″ north latitude). Yuikogane has a bright yellowish endosperm, a large grain size and is early heading (Nakajo 2015). Shinanotsubuhime was bred in Nagano prefecture, in the central region of Japan (36°39′05″ north latitude), and has a yellowish endosperm, a semi-dwarf habit and is late heading. Shinanotsubuhime was crossed with Yuikogane to create an F1 that was self-pollinated to generate the F2 population.

The 382 plants of the F2 generation and 10 plants of each parental cultivar were grown in a field in Iwate, Japan from June to October 2015, and the number of days to heading (DTH) under natural-day (ND) conditions was measured. The mean day-lengths during the cultivation period were 14.9 h in June, 14.6 h in July, 13.7 h in August, 12.4 h in September and 11.1 h in October. The mean temperatures were 12°C in June, 17.7°C in July, 18.4°C in August, 13.1°C in September and 4.1°C in October. In addition, the parental cultivars were grown under long-day (LD; 14 h light; 20000 lux, 25°C for 14 h and 22°C for 10 h) and short-day conditions (SD; 10 h light; 20000 lux, 25°C for 10 h and 22°C 14 h) in a controlled growth cabinet (NK system LH-220SP, Japan). We recorded the number of DTH for each plant as the number of days after sowing prior to the appearance of the first panicle.

Re-sequencing of Yuikogane and bulked samples for QTLseq

Genomic DNA was isolated from the leaves of Yuikogane using a DNeasy Plant Mini Kit (Qiagen) and used to construct the library for Illumina sequencing. The constructed libraries of Yuikogane were subjected to 250-bp paired-end sequencing using the Illumina Miseq (Illumina, CA, USA). The short reads of Yuikogane were aligned to the foxtail millet reference genome of a cultivar Yugu1 (Bennetzen ). For sequencing the bulk samples, we selected F2 individuals having a Yuikogane-type heading date (Y-progeny), a Shinanotsubuhime-type heading date (S-progeny) or a late-type heading date (L-progeny). See Results for number of individuals for each bulk type. Genomic DNA extracted from individuals of each progeny type was combined to make three bulk DNAs: a Y-bulk, an S-bulk and an L-bulk. The libraries of bulked DNAs were subjected to 75-bp paired-end sequencing using an Illumina NextSeq 500. Using the method of Takagi , we carried out QTL-seq analysis using short read sequences from each bulk. To obtain high-quality reads, short reads in which more than 20% of the sequenced nucleotides had a phred quality score of <20 were excluded from the analysis. The cleaned reads were aligned to the reference genome of Yuikogane using BWA aligner (version 0.7.15). After aligning the short reads, the Coval software (version 1.2) was used for filtering out mismatched reads and for variant calling (Kosugi ). The SNP-index was defined as the ratio between the number of reads of Shinanotsubuhime SNP and the total number of reads corresponding to the SNP (Abe ). We obtained SNP-index values of the bulked DNAs and calculated the Δ(SNP-index) whereby the Δ(SNP-index) = (SNP-index of Y-bulk) − (SNP-index of S-bulk or L-bulk). A sliding window analysis was applied by averaging the Δ(SNP-index) values within a 1 Mb window size and a 10 kb step increment.

QTL analysis with insertion-deletion polymorphisms (Indel) and cleaved amplified polymorphic sequence (CAPS) markers

To develop Indel and CAPS markers in the region near 43–44 Mb on chromosome 2 of Setaria italica, we searched for polymorphisms between the two parental lines by aligning Illumina reads to the reference genome of S. italica with BWA software (Bennetzen ). Primers for the Indel and CAPS markers were designed with Primer 3 (http://bioinfo.ut.ee/primer3-0.4.0/primer3). DNA was extracted from F3 seeds of each F2 individual by a CTAB method (Murray and Thompson 1980). For polymerase chain reaction (PCR) of Indel and CAPS markers, we used a 10 μl reaction volume containing 2 μl template DNA (10 ng μl−1), 0.05 μl of each primer (100 μM), 2.9 μl H2O and 5 μl Quick Taq HS DyeMix (Toyobo). The PCR profiles included an initial denaturation step at 94°C for 2 min followed by 40 cycles of 94°C for 30 sec, 58°C for 30 sec, 68°C for 15 or 30 sec, and finally an extension at 68°C for 5 min. To score CAPS genotypes, the amplified products were digested overnight in a 20 μl reaction volume containing 10 μl amplified product, 2 μl buffer, 7.8 μl H2O and 0.2 μl restriction enzymes, BamH I, BciT130 I, Bcn I, Hae III, Xsp I (Takara Bio), Dde I or Taq I (New England Biolabs) at the optimum reaction temperature for each enzyme. The amplified products of the Indel markers and digested products of the CAPS markers were electrophoresed in 3% and 1.5% agarose gel to detect polymorphisms, respectively.

Linkage maps were constructed using CarthaGene (version 1.2.3) with a Haldane mapping function (de Givry ). Markers were assigned to linkage groups using the “group” command with an LOD = 3.0 and a map distance below 30 cM. QTL analysis was carried out via composite interval mapping methods (CIM) with R/qtl (version 1.40-8) (Broman ). The threshold value (α = 0.05) for declaring the presence of a QTL was estimated by a 1000 times permutation test.

Results

Differences in photoperiodic sensitivity between Yuikogane and Shinanotsubuhime

The foxtail millet cultivars Yuikogane and Shinanotsubuhime differ in heading date and in their responses to photoperiod. The DTH for Yuikogane was 80.7 days under ND conditions, a value that was about 11 days shorter than that for Shinanotsubuhime (91.6 days) (Fig. 1A, 1B). To examine the effect of day-length on heading date for each cultivar, we grew the two cultivars under SD and LD conditions. Under SD conditions, the DTH for Yuikogane (27.8 days) was 17.3 days shorter than that of Shinanotsubuhime (45.1 days), whereas under LD conditions, the DTH for Yuikogane (43.6 days) was 31.9 days shorter than that of Shinanotsubuhime (75.5 days) (Fig. 1B). The DTH for Yuikogane was shorter than that of Shinanotsubuhime under all three growing conditions. The difference in the DTH for Yuikogane between SD and LD conditions (15.8 days) was smaller than that in Shinanotsubuhime (30.4 days). These results suggest a difference in the photoperiodic response between Yuikogane and Shinanotsubuhime.

Fig. 1

Days to heading (DTH) for Yuikogane and Shinanotsubuhime under different day-length conditions and frequency distributions of DTH in F2 population. A: The phenotype of Yuikogane (left) and Shinanotsubuhime (right) 85 days after sowing under ND conditions. The scale bars are 10 cm. B: DTH for Yuikogane and Shinanotsubuhime under natural day-length (ND) conditions in a field, long-day (LD) and short-day conditions (SD) in a growth cabinet. Error bars represent the standard error (n ≥ 10). C: DTH was investigated under natural day-length conditions in an F2 population (n = 381) of a cross between Shinanotsubuhime and Yuikogane. The mean values for DTH were 80.7 days (Yuikogane) and 91.6 days (Shinanotsubuhime) as indicated by white and black arrowheads, respectively. Red boxes indicate Y-bulk, S-bulk and L-bulk respectively. These bulk samples applied to QTL-seq.

Distribution of DTH in F2 plants derived from a cross between Shinanotsubuhime and Yuikogane

Under field conditions, transgressive segregation in DTH toward late heading was observed in F2 population (Fig. 1C). This transgressive segregation could not be explained by the effect of a single gene, indicating that multiple QTLs must be involved in the transgressive segregation of DTH in the F2 population. To carry out QTL-seq, we defined 22 individuals having early heading (DTH = 80 and 81 days) as Yuikogane-type (Y-) progeny, 32 individuals having late heading (DTH = 91) as Shinanotsubuhime-type (S-) progeny and 33 individuals having extremely late heading (DTH > 102) as late (L-) progeny.

Re-sequencing and QTL-seq analysis of DTH in foxtail millet

To construct a reference sequence for Yuikogane, we performed whole-genome re-sequencing using an Illumina MiSeq (DDBJ: DRR092734 and DRR092735). After aligning the obtained sequence reads to the Yugu1 reference genome (Bennetzen ), nucleotides of Yugu1 were replaced with those of Yuikogane at all SNP positions (1,102,448 positions) between the two varieties.

Each DNA bulk was subjected to whole-genome resequencing using an Illumina NextSeq. We obtained 67.3 million, 90.0 million and 65.9 million sequence reads from the bulk DNAs of Y-progeny, S-progeny and L-progeny (DDBJ: DRR089343, DRR089342 and DRR089341), respectively. When these reads were aligned to the developed Yuikogane reference sequence, the average depth was >6.39x for all bulked DNA (Supplemental Table 1), a sufficient depth for QTL-seq analysis (Takagi ). For QTL-seq, we made two comparisons, the “Y-bulk vs S-bulk” and the “Y-bulk vs L-bulk”.

To identify the candidate regions controlling the difference in DTH between Yuikogane and Shinanotsubuhime, we performed QTL-seq analysis of the “Y-bulk vs S-bulk”. A total of 45,370 SNPs was identified between Yuikogane and Shinanotsubuhime, and the SNP-index was calculated for each SNP. SNP-index plots of Y-bulk and S-bulk, and a Δ(SNP-index) plot for nine chromosomes are shown in Supplemental Fig. 1. We found contrasting patterns in the SNP-index graph of the Y-bulk and S-bulk in the region between 38.2 and 39.6 Mb on chromosome 2 as shown in Fig. 2A. The chance that the Δ(SNP-index) is lower than −0.49 as observed for the region of 38.2–39.6 Mb is P < 0.05 under the null hypothesis. Examining the SNP haplotype among early heading individuals in the Y-bulk showed that these plants carried the Yuikogane allele in the candidate region on chromosome 2, whereas late heading individuals in the S-bulk had the Shinanotsubuhime allele. These results indicated that there was a major QTL, named qDTH2, that controlled heading date within the 38.2–39.6 Mb region on chromosome 2. The Yuikogane allele of qDTH2 caused early heading and the Shinanotsubuhime allele caused late heading in the F2 population.

Fig. 2

QTL-seq applied to F2 population of a cross between Shinanotsubuhime and Yuikogane identifies QTLs for regulating DTH. The SNP-index was calculated based on 1 Mb interval with a 10 kb sliding window analysis (red line). Statistical confidence intervals for the null hypothesis of no QTLs (P < 0.05; green line). Red dotted boxes indicated candidate region identified by QTL-seq analysis. A: Genome-wide comparison of SNPs between the Y-bulk (early heading) and the S-bulk (late heading). SNP-index plots of the Y-bulk (top), S-bulk (middle) and Δ(SNP-index) (bottom) plots of chromosome 2. B: Genome-wide comparison of SNPs between the Y-bulk and the L-bulk. SNP-index plots of the Y-bulk (top), the L-bulk (middle) and the Δ(SNP-index) plot (bottom) of chromosome 7. C: Genetic linkage analysis with CAPS and Indel markers confirmed the location of qDTH2. D: Genetic linkage analysis with CAPS and Indel markers confirmed the location of qDTH7. Scale of y-axis shows lod value and scale of x-axis shows centimorgan (cM).

We performed QTL-seq analysis in the “Y-bulk vs L-bulk” to identify the candidate regions controlling transgressive segregation of extremely late heading in the F2 population. SNP-index plots of the Y-bulk and the L-bulk, and the Δ(SNP-index) plot for nine chromosomes are shown in Supplemental Fig. 2. We found contrasting patterns in the SNP-index graph for the Y-bulk and the L-bulk in the region between 29.2 and 31.0 Mb on chromosome 7 as shown in Fig. 2B. The chance that the Δ(SNP-index) becomes higher than 0.48 as observed for the region within 29.2–31.0 Mb is P < 0.05 under the null hypothesis. Observation of SNP haplotypes among the extremely late heading individuals in the L-bulk showed that most of these plants carried the Yuikogane allele in the candidate region on chromosome 7, whereas early heading individuals in the Y-bulk had the Shinanotsubuhime allele. These results suggested that there was a major QTL, named qDTH7, controlling extremely late heading within the 29.2–31.0 Mb region on chromosome 7.

Validation of the candidate QTLs detected by the QTL-seq method

To verify the candidate QTLs detected by QTL-seq analysis, we carried out QTL analysis using the CIM method for the same F2 population. We developed 24 markers, both Indel and CAPS markers, that were polymorphic between Yuikogane and Shinanotsubuhime in the candidate regions of chromosome 2 (36.3–41.8 Mb) and chromosome 7 (28.5–31.5 Mb) and a non-candidate region of chromosome 3 (46.5–50.5 Mb), respectively (Supplemental Table 2). Using these markers, we constructed a genetic linkage map for 214 F2 plants (Supplemental Fig. 3). A QTL with an LOD score = 3.49 was detected near Indel2_3 on chromosome 2 and another QTL with an LOD score = 13.96 was detected near CAPS7_3 on chromosome 7 (Fig. 2C, 2D); in contrast, no QTLs were detected on chromosome 3 (data not shown). The positions of these two QTLs corresponded to the genomic regions for qDTH2 and qDTH7 detected by the QTL-seq method. The phenotypic effect of qDTH2 was relatively small; the additive effect for the F2 population was 1.90 and the R value was 0.032 with the Shinanotsubuhime allele showing an increased effect on DTH. The phenotypic effect of qDTH7 was relatively large; the additive effect for the F2 population was −5.08 and R value was 0.486 with the Yuikogane allele showing an increased effect on DTH. These results were consistent with the results of the QTL-seq analysis.

We investigated the effect of different alleles of qDTH2 and qDTH7 between Yuikogane homozygous, heterozygous and Shinanotsubuhime homozygous plants (Fig. 3A, 3B). F2 plants homozygous for the Yuikogane allele at qDTH2 (Indel2_3) headed earlier (56.4 ± 0.9 days) than did those homozygous for the Shinanotsubuhime allele (60.0 ± 1.0 days). The number of DTH for heterozygous plants was intermediate between those of the homozygous plants (Fig. 3A). These results suggest that the Yuikogane allele at qDTH2 decreases the number of DTH in a semi-dominant manner and the effect of qDTH2 was approximately 4 days. F2 plants homozygous for the Yuikogane allele at qDTH7 (CAPS7_3) headed later (65.5 ± 1.0 days) than did those homozygous for the Shinanotsubuhime allele (55.4 ± 0.5 days). Heterozygous plants headed in a timeframe comparable (55.2 ± 0.2 days) to those homozygous for the Shinanotsubuhime allele, suggesting that the Yuikogane allele at qDTH7 increased the DTH in a recessive manner (Fig. 3B). Furthermore, to test whether genetic epistasis existed between qDTH2 and qDTH7, we compared the DTH among four genotype classes of the nearest marker (Indel2_3, CAPS7_3). A two-way ANOVA revealed the genetic effect of both QTLs in all genotype classes, although the interaction was not significant (p value = 0.77). These results suggest that the phenotypic variation in heading date was independently controlled by qDTH2 and qDTH7 in the F2 population.

Fig. 3

Validation of allelic effects of qDTH2 and qDTH7. A: The allelic effect of qDTH2 based on the genotypic classes at the Indel2_3 marker using 214 F2 individuals. B: The allelic effect of qDTH7 based on the genotypic classes at the CAPS2_3 marker using 214 F2 individuals. The asterisks indicate significant differences indicated by a Tukey-Kramer analysis. *P < 0.05, **P < 0.01.

Search for candidate genes in qDTH2 and qDTH7

To identify potential candidate genes in qDTH2 and qDTH7, the SNP index was calculated for all bulk samples and SNPs causing non-synonymous substitutions between the parents were selected. As qDTH2 might be a semi-dominant gene, we hypothesized that the range for the SNP index in the Y-bulk was 0–0.5 and the range for the S-bulk was 0.5–1. We found 68 SNPs to satisfy these requirements. Of these 68 SNPs, we found nine non-synonymous SNP and two deletions causing frame shifts in seven genes (Table 1). In qDTH7, SNPs were detected in 449 positions; SNPs of the L-bulk had a SNP index = 0 and the SNPs of the Y-bulk had a SNP index > 0.6 with a read depth of >5. Of the 449 SNPs, we found 44 SNPs causing a non-synonymous SNP, one SNP altering a start codon, and single deletions causing a frame shift in 31 genes (Table 2). Although there have been two genes reported to be involved in heading and flowering in rice, maize, and Arabidopsis within the corresponding candidate regions of qDTH2, we did not detect SNPs causing non-synonymous substitutions in these genes. In the candidate region for qDTH7, we found that Seita.7G246700, a homolog of Rice outermost cell-specific gene (Roc4), had non-synonymous substitutions by two SNPs. Seita.7G246700 is predicted to encode a homeobox domain and a START domain. Comparison of the amino acid sequences of Seita.7G246700 from Yuikogane and Shinanotsubuhime revealed a non-synonymous substitution in a START domain.

Table 1

Identification of SNPs in putative candidate genes around the qDTH2

Gene name	Positions (bp)	Reference	Variant	Variant effect	Y-bulk variant rate	S-bulk variant rate	Description
Seita.2G285800	38231339	–	A	frame shift	0.14	1	No protein domain
Seita.2G290100	38580942	G	A	missense	0.17	0.69	Protein of unknown function
Seita.2G293600	38903157	A	G	missense	0.4	0.82	Transferase family
Seita.2G296100	39049199	–	GC	frame shift	0.25	0.88	GLUTATHIONE S-TRANSFERASE
Seita.2G296300	39075668	A	C	missense	0.4	0.6	AUXIN-RESPONSIVE FAMILY PROTEIN
Seita.2G296500	39084573	G	A	missense	0	0.71	AUXIN-RESPONSIVE FAMILY PROTEIN
	39084644	G	C	missense	0.4	0.73
Seita.2G297100	39107909	A	G	missense	0	1	REGULATOR OF VPS4 ACTIVITY IN THE
	39108137	C	T	missense	0.33	1	MVB PATHWAY PROTEIN
	39108296	C	T	missense	0.25	1
	39109931	C	A	missense	0	0.67

Table 2

Identification of SNPs in putative candidate genes around the qDTH7

Gene name	Positions (bp)	Reference	Variant	Variant effect	Y-bulk variant rate	L-bulk variant rate	Description
Seita.7G232300	29604317	C	G	missense	1	0	OPT oligopeptide transporter
Seita.7G232400	29606996	A	G	missense	0.7	0	OPT oligopeptide transporter
	29607005	G	A	missense	0.75	0
Seita.7G232400	29607117	G	A	missense	0.66	0	AN1-TYPE ZINC FINGER PROTEIN
Seita.7G236700	29889610	G	A	missense	0.88	0	COPPER TRANSPORT PROTEIN ATOX1
Seita.7G236900	29903348	G	T	missense	0.33	0	LEUCINE RICH REPEAT
	29903684	T	G	missense	0.71	0
	29903698	C	T	missense	0.71	0
	29904153	G	T	missense	1	0
	29904176	C	T	missense	1	0
	29904227	C	T	missense	0.83	0
	29905361	C	G	missense	0.8	0
Seita.7G237100	29913322	C	T	missense	0.86	0	No domain
Seita.7G237300	29927222	A	C	missense	0.86	0	Protein tyrosine kinase
Seita.7G237600	29974234	C	G	missense	0.78	0	Non-specific protein-tyrosine kinase
Seita.7G237800	30013535	AG	GC	missense	0.83	0	Protein kinase domain (Pkinase)
	30041257	G	A	missense	0.88	0
	30041336	G	A	missense	0.8	0
Seita.7G238500	30070539	G	A	missense	0.7	0	SEED STORAGE 2S ALBUMIN SUPERFAMILY PROTEIN
Seita.7G238600	30074953	A	G	missense	0.67	0	PEROXISOME ASSEMBLY FACTOR
	30076086	A	C	missense	0.67	0
Seita.7G238800	30093672	T	C	missense	0.88	0	IMIDAZOLEGLYCEROL PHOSPHATE DEHYDRATASE HIS7
Seita.7G239000	30102013	A	T	missense	0.71	0	PPR repeat family (PPR_2)
Seita.7G239300	30116522	G	C	missense	0.71	0	No domain
Seita.7G239400	30132361	G	C	missense	0.73	0	No domain
Seita.7G239500	30134969	–	deletion	frame shift	0.75	0	F-box domain (F-box)
Seita.7G239600	30139521	T	C	start lost	0.71	0	No protein domain
Seita.7G240300	30191072	A	C	missense	1	0	Kelch motif
	30297454	T	C	missense	0.67	0
Seita.7G241700	30306287	C	G	missense	0.89	0	Protein of unknown function
Seita.7G243900	30427428	G	C	missense	0.67	0	GLYCOSYLTRANSFERASE
Seita.7G244100	30438355	G	C	missense	0.75	0	Cysteamine dioxygenase/Persulfurase
Seita.7G244300	30447398	C	T	missense	0.67	0	F11F12.2 PROTEIN
	30448527	C	A	missense	0.6	0
Seita.7G245700	30512205	A	T	missense	0.8	0	No domain
Seita.7G245900	30519724	C	G	missense	1	0	LEUCINE-RICH REPEAT-CONTAINING PROTEIN
Seita.7G246700	30568896	A	C	missense	0.7	0	HOMEOBOX-LEUCINE ZIPPER PROTEIN HDG2
	30569082	T	C	missense	1	0
Seita.7G249600	30709399	T	G	missense	0.71	0	Anthocyanidin reductase
Seita.7G249700	30717256	G	A	missense	0.78	0	Anthocyanidin reductase
	30718052	G	C	missense	0.75	0
Seita.7G249800	30725051	C	G	missense	0.75	0	Anthocyanidin reductase
Seita.7G249900	30730621	T	C	missense	0.75	0	ORGANIC CATION/CARNITINE TRANSPORTER 4
Seita.7G250000	30739415	–	deletion	missense	0.71	0	ATP-BINDING TRANSPORT PROTEINRELATED
Seita.7G250400	30774875	C	G	missense	0.67	0	Protein of unknown function
	30774997	AA	CC	missense	0.75	0

Discussion

In this study, we used an F2 population from a cross between Yuikogane and Shinanotsubuhime and identified two QTLs, qDTH2 and qDTH7, that are associated with heading time in foxtail millet. In the conventional approaches, construction of a genetic linkage map and QTL analysis have been required to develop molecular markers and genotyping every individual in a mapping population. On the other hand, the QTL-seq method applied here only required whole genome re-sequencing of DNAs from two or more bulked samples that have extremely different traits in the segregating progeny, reducing cost and efforts. As parents of the F2 population were both inbred lines from Japanese cultivars; polymorphisms between these parents were sufficient to be detected by QTL-seq analysis. Photoperiodic sensitivities were different between Yuikogane and Shinanotsubuhime (Fig. 1B), and QTL qDTH2 was identified by QTL-seq comparison of bulked samples, a Y-bulk and an S-bulk (Fig. 2). These results indicate that qDTH2 regulates both differences in photoperiodic sensitivity and DTH in Yuikogane and Shinanotsubuhime. Also, qDTH2 may be related to the natural variation of DTH among Japanese cultivars.

A two-way ANOVA revealed that there is no epistasis between qDTH2 and qDTH7; however, extremely late heading was observed in only 33 individuals with the Yuikogane allele of qDTH7 in F2 segregating progeny (n = 381). In general, a transgressive segregation pattern is caused by allelic interaction between parents. Therefore, we hypothesize that there is at least one other Shinanotsubuhime allele involved in the extremely late heading phenotype that was not detected by QTL-seq analysis, because this additional allele may function as dominant. Further analysis is needed to identify the gene(s) responsible for the extremely late heading using RILs from Yuikogane and Shinanotsubuhime.

As a few research groups have reported that the homolog of rice HD1 gene is a candidate of a QTL for DTH by QTL analysis and GWAS in foxtail millet (Jia , Mauro-Herrera , Zhang ), we compared HD1 gene sequence between Yuikogane and Shinanotsubuhime based on the aligned short reads of NGS. However, we did not found any polymorphism between two parental cultivars. Liu have performed DNA sequence analysis of HD1 orthologs from 15 wild and 83 domesticated accessions in foxtail millet, and found 1 splice site from “GT” to “AT” in first exon, “GT” type and “AT” type are high assosiation with wild and domesticated accessions, respectatively. Yuikogane and Shinanotsubuhime have HD1 allele of domesticated “AT” type. Taken together, these findings suggest that HD1 gene does not contribute to variation of DTH in the F2 population.

The position of qDTH2 was defined in a specific genomic region (38.2–39.6 Mb) by QTL-seq and near the marker Indel2-3 (38.4 Mb) by using the CIM method. qDTH2 is located near QTL 2.2 that was previously reported by Mauro-Herrera to be a QTL controlling heading date. QTL2.2 located in the genomic region of 35.5–39.0 Mb on foxtail millet chromosome 2 was defined by two flanking markers, UGSF242 and UGSF249, and was closest to UGSF248 (38.3 Mb). This region contains six candidate homolog genes associated with DTH in rice and maize, whereas qDTH2 contains two common genes (Seita.2G286100 and Seita.2G291300) that are homologous to Oryza sativa Pseudo-Response Regulator95 (OsPRR95) and Zea mays Delayed flowering1 (DLF1). OsPRR95 encodes a highly homologous PRR protein and plays a role in the circadian clock with OsPRR1, OsPRR37, OsPRR59 and OsPRR73 in rice (Murakami ). Of these PRR genes, OsPRR37 functions as a major long-day dependent flowering repressor with grain number, plant height, and heading date7 (Ghd7) and plays an important role in photoperiod sensitivity in rice (Kim , Koo , Kwon ). DLF1 is similar to FLOWERING LOCUS D (FD) in Arabidopsis; both genes encode a basic leucine zipper protein expressed in the shoot apex to induce the floral transition (Abe , Muszynski ). However, the two genes (Seita.2G286100 and Seita.2G291300) from Yuikogane and Shinanotsubuhime are identical in their coding regions. Among the 139 genes in the genomic region of qDTH2, we found non-synonymous SNPs in seven genes (Table 1), but none of these genes seem to be involved in heading based on their description and sequence homology with Arabidopsis, rice and maize genes. Because we extracted SNPs and small Indels from the Illumina short reads of genomic DNA, it is not possible to detect relatively large Indels or differences in transcription products between the parents. Because we focused on only non-synonymous SNPs, we may have overlooked mutations in other candidate genes or there is a possibility that causal genes are located outside of the region. Further evidence is needed to identify genes regulating heading in the candidate region of qDTH2.

The position of qDTH7 was defined in a specific genomic region (29.2–31.0 Mb) by QTL-seq and near the marker CAPS7_3 (30.4 Mb) based on the CIM analysis. The qDTH7 is located near QTL7.1, another QTL reported by Mauro-Herrera that controls heading date. QTL7.1 is in the 31.3–34.0 Mb region of foxtail millet chromosome 7 and is flanked by the two markers, UGSF665 and UGSF778. Seita.7G263000, a homolog of APETLA 2 (AP2), resides near qDTH7 and QTL7.1, and there are non-synonymous SNPs between Yuikogane and Shinanotsubuhime alleles. However, these are not located in the AP2 domain. The location of Seita.7G263000 gene is 1.1 Mb away from the peak position of qDTH7 (CAPS7_3), suggesting that this gene is not responsible for the extremely late heading date and that qDTH7 and QTL7.1 are close but different QTLs.

Among the 257 genes predicted in the qDTH7 region, we found non-synonymous substitutions by two SNPs in the Roc4 homolog, Seita.7G246700. Roc4 promotes flowering time under long days in rice (Wei ), suggesting that the foxtail millet homolog of Roc4 may be a possible candidate for qDTH7. Roc4 is similar to OUTER CELL LAYER1 (ZmOCL1) in maize and protodermal factor 2 (AtPDF2) in Arabidopsis, both of which encode a homeodomain leucine zipper (HD-zip) class IV family protein (Abe , Depege-Fargeix , Wei ). Several HD-zip IV genes are thought to be related to flowering time. AtPDF2 regulates flowering, and its overexpression delays flowering (Abe ). Furthermore, ZmOCL1 suppresses the floral transition, but Roc4 activates flowering under long day conditions (Depege-Fargeix , Wei ). Roc4 RNAi plants, AtPDF2 null mutants and ZmOCL1 RNAi plants have been reported not to change their flowering phenotype compared with the wild-type plants (Abe , Depege-Fargeix , Wei ), suggesting that these are probably functionally redundant genes. Therefore, qDTH7 may also be a redundant QTL and regulates the extremely late heading phenotype together with other QTLs.

In conclusion, we detected two foxtail millet QTLs for heading date, qDTH2 and qDTH7, using QTL-seq. We confirmed the QTLs with conventional linkage analysis in the candidate region using CAPS and Indel markers that were developed from the genomic sequence obtained by NGS. The allelic difference in DTH at qDTH2 was relatively small (about 4 days). Such a QTL with a small effect is very important for modifying foxtail millet DTH since slight changes in flowering time could be of value to breeders. We need to experiment further with QTL fine mapping of these newly identified QTLs using RILs or near isogenic lines (NILs). Furthermore, RNA-seq analysis will facilitate our understanding of the transcription network responsible for variation in DTH in foxtail millet.