The seeds of foxtail millet variety “Yu Gu 18” were planted and grown in the experimental station of Anyang Institute of Technology (Anyang, Henan, China). The same-aged panicle with the same size was harvested from T1–T5 developmental stages (T1: 7 days after anthesis, T2: 14 days after anthesis, T3: 21 days after anthesis, T4: 28 days after anthesis, and T5: 35 days after anthesis). The experimental field, agronomical practices situation, and sampling process are described in the previous study (Wang et al., 2020).

In isolating genomic DNA, DNeasy Blood and Tissue Kit (Qiagen, Dusseldorf, Germany) was used. Fifteen high-quality BS libraries were constructed and then sequenced with the Illumina Hiseq X Ten platform by Shanghai OE Biotech (Shanghai, China). Briefly, DNA was first fragmented to an approximate size of 250 bp. Then, dA was added to the 3′ end by blunt-end cloning and methylated adaptor ligation. The EZ DNA Methylation-Gold kit (ZYMO, Tustin, CA, United States) was used to bisulfite-convert the ligated DNA. Different DNA fragment lengths were excised from a 2% agarose gel and purified and amplified by PCR.

All sequence data were deposited into the NCBI database (Accession number: PRJNA699635). Clean data were obtained by filtering the low-quality sequences and mapped to the foxtail millet genome (Bennetzen et al., 2012) using the mapping program, Bismark (Krueger and Andrews, 2011).

Differentially methylated regions (DMRs) were identified by MethylKit (Akalin et al., 2012), which used a sliding window approach. The window was 1,000 bp, and the step length was 100 bp. The logistic regression was applied to detect significant DMRs and the screened criteria as methylation difference of ≥10% and p ≤ 0.05. After DMRs were identified, differentially methylated genes (DMGs) located in DMRs were characterized.

For the identification of differentially methylated promoter (DMP), 2 kb upstream of the transcription start site (TSS) was used as the gene promoter region, and the methylation level in this region was counted. The methylation level of the promoter region between different samples was compared. To analyze DMP between groups, a t-test was used. The differences in methylation levels of all genes among all groups were analyzed. DMP regions were screened according to methylation differences of ≥10% and p ≤ 0.05.

The total RNA of kernels from spikelets was extracted using a mirVana miRNA Isolation Kit (Ambion, Inc., Austin, TX, United States) with three biological replicates. Afterward, high-quality libraries were sequenced on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, United States) with 125 paired-end sequencings. The data analysis was performed according to previously described methods (Wang et al., 2020).

For integration with DNA methylation, RNA-seq data were used representing the same stages of grain filling in foxtail millet from the previous study.

The Pearson correlation analysis of RNA-seq data and DNA methylation levels in each sample was completed as previously described (Jin et al., 2014).

Differentially expressed genes associated with DMR were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis¹ (Kanehisa et al., 2008) to map them to the specific pathways they were involved in using the Blast_v2.2.26 software.

Deoxyribonucleic acid libraries for five different stages were constructed for sequencing to get the characteristics and change of DNA methylation during grain development. A total of 2.16 billion raw reads were obtained, and the mapping rate was 76.6–83.5%. Moreover, conversion rates of each library were >99.56% (Supplementary Table 1). Sequencing depth was >30-fold coverage per DNA strand (Supplementary Table 2). Furthermore, the proportion of the 5 mC methylation for three contexts in each chromosome was shown in Supplementary Figure 1. The genome-wide pattern of methylated DNA was similar across chromosomes at different stages (Figure 1A and Supplementary Figure 2). The methylation level of different chromosomal positions showed significant differences, and the lowest CHH methylation was found across the genome (Figure 1B and Supplementary Figure 3).

The context CHH methylation was accounted for the largest percentage during the whole grain filling stage followed by CG and CHG methylation. Furthermore, CHH methylation remained elevated during grain filling (38–46%), but CG and CHG methylation showed the opposite trend (33–29 and 29–25%, respectively) (Figure 2A). The mean methylation levels were lowest at T2, followed by T4, which mainly occurred in the CHG context (Figure 2B). Moreover, methylation level in the promoter, gene-body, and downstream regions of genes also reached the bottom at T2, and then at T4 (Figure 2C). Besides, methylation patterns of TEs were compared among different stages. No significant differences were found in the 2 kb upstream region among five developmental stages. Interestingly, the CG methylation level of the TE-body region was less at T2 than that at other stages (Figure 2C). The CHH methylation level of the TE-body and 2-kb downstream region was methylated the least at T1 (Figure 2C). Circos plots of DNA methylation at later stages compared with T1 uncovered the distinctions in the three contexts at different genomic regions (Figure 3 and Supplementary Figure 4).

The differentially methylated region was analyzed for understanding the methylation differences at later stages compared with T1. The number of hyper-DMRs identified at later stages relative to T1 was gradually increased from T2 to T3, then decreased from T3 to T4, and increased again from T4 to T5. The hypo-DMRs showed an opposite trend. Furthermore, the number of CG DMRs was far more than other types of DMRs (Supplementary Figure 5). In CG and CHG context, the exon and intergenic regions had the most DMRs, respectively. In the CHH context, the promoter regions had the most hypo-DMRs, while hyper-DMRs were mostly distributed in the intergenic (T2 and T4), exon (T3), and promoter (T5) (Supplementary Figure 6). The numbers of genes with DMRs in gene-body regions (DMR_genes) and promoter regions (DMP_genes) among the different comparisons were shown in Venn analysis. For example, a total of 7,460, 2,186, and 907 DMR genes were identified for CG, CHG, and CHH contexts in the T2 vs. T1, and the number of DMP genes was 1,234, 1,234, and 183 (Figure 4A). Additionally, the methylation levels of the DMRs were decreased at later stages, but T3 and T5 showed higher CHH methylation levels than those in the T1 stage (Figure 4B).

A total of 11 methyltransferases, including 3 chromomethylases (CMT), 4 domain-rearranged methyltransferases (DRM), 3 methyltransferase 1 (MET1), and 1 DNA methyltransferase-2 (DNMT2), were identified (Supplementary Table 3). Among the 11 methyltransferase genes, CMT3, DRM3, MET1-1, and MET1-3 are barely expressed (FPKM < 1) at all five stages (Supplementary Table 3). The foxtail millet genome harbors 7 DNA demethylase homologs, including 2 DEMETER (DME), 1 DEMETER-like2 (DML2), DEMETER-like3 (DML3), and 3 Repressor of Silencing1 (ROS1). DML2 and ROS1-3 are barely expressed (FPKM < 1) at each developmental stage (Supplementary Table 3).

A correlation analysis was performed to analyze the relationship among DNA methylation levels, DNA methyltransferase, and demethylase (Supplementary Table 4). The expression level of SiCMT2 and SiDRM4 were highly positively correlated with methylation levels of C, CG, and CHG, but showed a highly negative correlation with CHH. The expression level of SiDRM2 was highly positively correlated with methylation levels of CG but showed a highly negative correlation with C, CHG, and CHH. The expression level of SiMET1-2 was highly positively correlated with methylation levels of C and CHH.

The deoxyribonucleic acid demethylase gene SiDME1 showed a highly positive correlation between expression and methylation level of C and CG. The expression level of SiDME2 was highly positively correlated with methylation levels of C and CHG but showed a highly negative correlation with methylation levels of CHH. The expression levels of SiDML3 and SiROS1-1 were highly positively correlated with the methylation levels of CHH and C.

Based on the FPKM, all genes were classified into four groups (Figure 5). Genes in the none-group were highest CG methylated in all genetic regions (Figure 5). The levels of none-expression genes of CHG methylation were highest in gene-body and upstream. By contrast, genes in the none-group were the lowest CHH methylated in all genetic regions. As expected, high expression genes had the lowest CG and CHG methylation and the highest levels of CHH methylation in gene-body. In low and moderate groups, no significant correlation was found between gene expression and methylation status of CG and CHH in most genetic regions. By contrast, genes with higher expression levels presented lower CHG methylation levels in gene-body (Figure 5).

Methylated genes were also grouped by the promoter and gene-body methylation levels (Figure 6). The methylation level in the gene-body region was positively correlated with the gene expression level, except the high methylated genes (fifth group), which had the lower expression levels (Figure 6A). However, the methylation level in the promoter region of most genes had no significant effect on gene expression, except the genes in the first group that also showed a lower expression level (Figure 6B).

To verify the effect of altered DNA methylation on DEGs, we focused on the overlapped genes of DMGs and DEGs for different comparisons. The correlation between DEGs and methylation was calculated by Pearson correlation analysis based on the expression level of DEGs and methylation level of DMR. We identified the DEGs associated with DMR in line with Pearson correlation ≥ 0.90 (positive correlation) or ≤–0.90 (negative correlation). There were 77, 49, and 14 DEGs associated with CG, CHG, and CHH DMR for T2 vs. T1, respectively (Figure 7B). DEGs associated with CG DMR were significantly enriched in photosynthesis-antenna proteins, metabolic pathways, and cysteine and methionine metabolism, respectively (Figure 7A). Representative CG DMRs in the promoter of LHCB5 (which encoded chlorophyll a-b binding protein CP26) were methylated less at T2 than at T1 and expressed at significantly lower levels at T2 than at T1.

The enriched categories of DEGs associated with CHG DMR included amino sugar and nucleotide sugar metabolism, biosynthesis of secondary metabolites, and phenylalanine metabolism (Figure 7A). Representative DMRs were methylated less at T2 than that at T1, such as APS1 (which encoded glucose-1-phosphate adenylyltransferase small subunit) and VTC1 (which encoded mannose-1-phosphate guanylyltransferase 1) in the CHG context in the promoter, whereas an opposite trend was observed in the expression level (Figure 7C).

Differentially expressed genes associated with CHH DMR were enriched in fructose and mannose metabolism, inositol phosphate metabolism, and protein export. In addition, representative CHH DMRs of VTC1 were methylated less at T2 than at T1, also expressed at significantly higher levels at T2 than at T1.

For T3 vs. T1, 245, 137, and 24 DEGs associated with CG, CHG, and CHH DMR were identified, respectively (Supplementary Figure 7B). KEGG enrichment of DEGs associated with CG DMR showed that the four most significantly altered pathways were metabolic pathways, glyoxylate, and dicarboxylate metabolism, inositol phosphate metabolism, and starch and sucrose metabolism, respectively (Supplementary Figure 7A). Representative CG DMRs in promoter or gene-body of four DEGs associated with starch and sucrose metabolism, SS4 (which encoded soluble starch synthase), PHO (which encoded alpha-1,4 glucan phosphorylase L isozyme), DPE1 (which encoded 4-alpha-glucanotransferase), and BGLU31 (which encoded beta-glucosidase 31) were methylated less at T3 than at T1 but expressed at significantly higher levels at T3 than at T1. In contrast, the expression level of FRK6 (which encoded fructokinase-6) was positively correlated with its methylation (Supplementary Figure 7C).

Differentially expressed genes associated with CHG DMR were enriched in seleno compound metabolism, other glycan degradation, and biosynthesis of secondary metabolites, respectively (Supplementary Figure 7A). Representative CHG DMRs in the promoter of MS1 (which encoded Cobalamin-independent methionine synthase 1) were methylated less at T3 than that at T1, and also expressed at significantly lower levels at T3 than at T1.

Differentially expressed genes associated with CHH DMR were enriched in porphyrin and chlorophyll metabolism, MAPK signaling pathway, and plant hormone signal transduction (Supplementary Figure 7A). The methylation of CHH DMR in gene-body of CHLH (which encoded magnesium-chelatase subunit) showed a negative correlation with its expression level.

For T4 vs. T1, 150, 134, and 34 DEGs associated with CG, CHG, and CHH DMR were identified respectively (Supplementary Figure 8B). KEGG pathway analysis of DEGs associated with CG DMR indicated that the three most significantly changed pathways were pyrimidine metabolism, oxidative phosphorylation, and purine metabolism (Supplementary Figure 8A). Representative CG DMRs in gene-body of four DEGs associated with pyrimidine metabolism, POLA (DNA polymerase I A), pyrH (uridylate kinase), RNR1 (ribonucleoside-diphosphate reductase large subunit), and PBY1 (5′-nucleotidase SurE) were methylated less at T4 than at T1, and the majority of them were expressed at significantly lower levels at T4 than at T1 except for POLA (Supplementary Figure 8C).

Kyoto Encyclopedia of Genes and Genomes pathway analysis demonstrated that the DEGs associated with CHG DMR were primarily related to photosynthesis, beta-alanine metabolism and ascorbate, and aldarate metabolism (Supplementary Figure 8A). Representative CHG DMRs in promoter and gene-body of FDX5 (which encoded ferredoxin) and FDC1 (which encoded ferredoxin-2) were methylated less at T4 than at T1 and were expressed at significantly lower levels at T4 than at T1.

Differentially expressed genes associated with CHH DMR were enriched in the biosynthesis of unsaturated fatty acids, alpha-linolenic acid metabolism, and pentose phosphate pathway (Supplementary Figure 8A). The methylation of CHH DMR in the promoter of ACX4 (which encoded acyl-coenzyme A oxidase 4) showed a positive correlation with its expression level.

There were 249, 147, and 28 DEGs associated with CG, CHG, and CHH DMR identified at T5 compared with T1, respectively (Supplementary Figure 9B). The top three KEGG pathways of DEGs associated with CG DMR were nitrogen metabolism, base excision repair, and anthocyanin biosynthesis (Supplementary Figure 9A). Representative CG DMRs in promoter or gene-body of DEGs associated with nitrogen metabolism, ACA7 (alpha carbonic anhydrase 7) and GLN2 (glutamine synthetase) were methylated and expressed less at T5 than that at T1, but the expression level of NRT2.4 (high-affinity nitrate transporter 2.4) was negatively correlated with its methylation.

The functional enrichment of DEGs associated with CHG DMR primarily included arginine and proline metabolism, phenylalanine, tyrosine, and tryptophan biosynthesis, and tropane, piperidine, and pyridine alkaloid biosynthesis (Supplementary Figure 9A). Representative CHG DMRs in the promoter of PAO1 (which encoded polyamine oxidase 5) and ASP5 (which encoded aspartate aminotransferase) were methylated less at T5 than that at T1 and expressed at significantly lower levels at T5 than at T1.

Differentially expressed genes associated with CHH DMR were enriched in the pentose phosphate pathway, fructose, and mannose metabolism, and MAPK signaling pathway (Supplementary Figure 9A). The methylation levels of CHH DMR in gene-body and promoter of PFP-ALPHA (which encoded pyrophosphate–fructose 6-phosphate 1-phosphotransferase subunit alpha) and G6PDH (glucose-6-phosphate 1-dehydrogenase) were higher at T5 than at T1 but expressed at significantly lower levels at T5 than at T1.