Comparative transcriptome analysis provides insights into grain filling commonalities and differences between foxtail millet [Setaria italica (L.) P. Beauv.] varieties with different panicle types

Plant materials and sample collection

The varieties used in our study were “Yu Gu 18” and “Zhang Gu 13”, grown in the experimental station of Anyang Academy of Agriculture Sciences (North China summer-sowing region, N 36°18′, E 114°37′, altitude 75 m) from May to October 2020. A randomized block design was used to perform the field experiments, with 20 m² for each plot and each plot receiving an appropriate supply of fertilizer and water during the growth period. Fertilizer was applied in the following amounts: N: 544 kg/hm², P₂O₅: 489 kg/hm², K₂O₅: 77 kg/hm². Starting 7 days after anthesis, same-aged panicles of similar sizes were harvested every 7 days, up to day 42 (T1–T6). Five panicles were dried at 105 °C for 0.5 h and 70 °C until their weight was constant, respectively. Panicles were weighed at the T1, T2, T3, T4, T5, and T6 stages of development. The significance among the means of panicle weight was verified using ANOVA, and mean differences were compared by Fischer’s least significant difference (LSD) test (p < 0.05).

RNA extraction, library construction, and sequencing

The total RNA of the kernels of each panicle was extracted using TRIzol Reagent (Ambion, Inc., Austin, TX, USA) according to the manufacturer’s protocol. RNA purity and yield were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The libraries were then constructed using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. The RNA was sent to Shanghai OE Biotech for library construction and 125 paired-end sequencing using an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA). All sequence data were deposited into the NCBI database (Accession number: PRJNA721299, SAMN18713804 (YG), SAMN18713832 (ZG)).

Data quality control and read mapping

Deep sequencing data were processed using Trimmomatic version 0.36 (Bolger, Lohse & Usadel, 2014). Clean reads were obtained by removing the low-quality sequences. Hisat2 version 2.2.1.0 (Kim, Langmead & Salzberg, 2015) was subsequently used to map the clean reads to the foxtail millet reference genome version 2.2 (Setaria italica v2.2: https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Sitalica).

Temporal expression patterns of genes

We employed the Short Time-series Expression Miner (STEM) software (Carnegie Mellon University, Pittsburgh, PA, USA) to analyze the temporal expression patterns of genes during grain filling. The gene expression data of each sample was imported (three replications), and the default settings of the STEM program were used. A hyper geometric distribution test was also carried out to identify the significantly enriched profiles (p value < 0.05) (Ernst & Bar-Joseph, 2006).

Differential expression analysis of genes

Principal component analysis (PCA) was performed to assess variability between the two assessed accessions using the PRCOMP command with default setting in the R-software package. The gene expression values were calculated as fragments per kilobase per million reads (FPKM), using Cufflinks version 2.2.1 (Trapnell et al., 2010). Differentially expressed genes (DEGs) were identified using the DESeq version 1.18.0 software with the R package functions ‘estimateSizeFactors’ and ‘nbinomTest.’ A gene was deemed differentially expressed with the absolute value of fold change > 2 and a p-value < 0.05 (Anders & Huber, 2012). The Venn and UpSet diagrams were drawn using TBtools with default settings (Chen et al., 2020).

Functional annotation of genes

All DEGs were subjected to a Gene Ontology (GO, http://www.geneontology.org/) (Kanehisa et al., 2008) and functional enrichment analysis to assign their putative function. In addition, we performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (http://www.genome.jp/kegg/pathway.html) (Ashburner et al., 2000) to map DEGs to the specific pathways that each gene is specifically involved in using the Blast_v2.2.26 software.

Validation of RNA sequencing data

The expression profiles of 13 randomly selected genes were confirmed through qRT-PCR (Table S1). The RNA of the kernels was extracted using TRIzol reagent according to the manufacturer’s protocol (Tiangen, Beijing, China). The Primer-BLAST online NCBI tool was used to design primers specific to each DEG. The cDNAs were synthesized using a FastKing RT kit (Tiangen, Beijing, China) and qRT-PCR was performed following the protocol of the SYBR Green PCR Master Mix Reagent (SuperReal PreMix Plus, Tiangen, China) on an ABI 7500 fast Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). Seita.7G294000 (Actin 7) was used as the internal reference (Qin et al., 2020; Yu et al., 2020) and the 2^−∆∆Ct method was used to calculate the relative expression level of each assessed gene (Livak & Schmittgen, 2001).

Panicle weight variations during grain filling

Panicles were harvested at six different time points (T1–T6) to compare the changes in panicle weight between the two assessed foxtail millet accessions. Compared to the YG accession, the panicle weight of ZG plants was significantly higher during the early stages (T1 to T3) of reproductive development (Fig. 1A). However, the panicle weight of YG increased rapidly and exceeded that of ZG during the late stages (T4 to T6) of panicle development (Fig. 1A).

Transcriptome analysis of YG and ZG

Thirty-six libraries were constructed from the six assessed developmental stages for RNA-seq to obtain the gene expression dynamics of the two assessed foxtail millet accessions during grain filling. A total of 904 and 1,026 million raw reads were obtained for YG and ZG, respectively (Table S2). Subsequent removal of the adaptor sequences and low-quality sequences yielded 887 and 1,005 million clean reads with a Q30 base percentage higher than 94.89% and 94.38% for YG and ZG, respectively. Read mapping revealed 89.94% of YG reads mapped to the genome, with a unique mapping rate of 83.86%. Similarly, the mapping and unique mapping rates of ZG reached 91.81% and 88.08%, respectively.

The relative global relationships of biological replicates and clustering patterns of the samples were assessed through principal component analysis (PCA). Three samples of each time point clustered together (Fig. 1B). YG samples were separated from ZG samples, independent of the different developmental stages, suggesting that the genetic backgrounds are indeed different.

Normalized read counts for each gene was calculated and divided into four groups based on their expression level (Fig. 1C). Genes identified in all stages of the same accession were combined, and a Venn diagram was used to reveal the genes that were commonly expressed between the YG and ZG (Fig. 1D). Notably, there were 10,479 commonly expressed genes between the YG and ZG during the six assessed time points as part of the grain filling process. Moreover, there were 347 genes that were only expressed in YG, and 267 expressed genes that were only identified in ZG.

Gene expression analysis of YG and ZG

The STEM (Short Time-series Expression Miner) software was used to cluster the genes which were expressed across the grain filling process in the two assessed accessions to determine their temporal expression patterns. Twelve and fifteen significantly enriched profiles (p-value < 0.05) were identified in YG and ZG from fifty profiles, respectively. Genes with similar expression trends were grouped into the same profile. Two particular model profiles (profile 8, downregulated profile; profile 39, upregulated profile) were identified (Fig. 2A). Venn diagrams revealed that both profiles had numerous uniquely expressed genes (Fig. 2B), suggesting that each accession utilises a unique regulation network during the grain filling progress. A KEGG enrichment analysis of profile 8 and profile 39 genes in both accessions revealed that genes involved in photosynthesis were downregulated in YG across the grain filling process, while those involved in ribosome biogenesis and the spliceosome were upregulated. Figure 2C shows the top 10 enriched pathways. In ZG, genes involved in starch and sucrose metabolism, phenylpropanoid biosynthesis and photosynthesis were downregulated, while those involved in peroxisome function and zeatin biosynthesis were upregulated across the grain filling process.

Pairwise comparisons were performed to investigate the transcriptome dynamics during grain filling in YG. We identified 6,388, 1,929, 922, 1,056 and 886 DEGs in different YG comparison sets (Fig. 3A). A Gene Ontology (GO) enrichment analysis of the DEGs was performed to analyze the DEGs based on their biological process (BP), cellular component (CC), and molecular function (MF). Figure S1 shows the top 10 BP, CC, and MF GO enrichment terms of each comparison. Figure 3B shows the significantly enriched pathways among the different comparison sets. Similarly, there were 4,485, 3,259, 1,058, 1,451 and 2,366 DEGs identified in different comparison sets of ZG (Fig. 3A). Figure S2 shows the top 10 BP, CC, and MF GO enrichment terms of the DEGs during grain filling in ZG. Figure 3B shows the significantly enriched pathways among different ZG comparison sets.

Identification of DEGs between YG and ZG

We identified 6,832 unique DEGs between YG and ZG (Data Set S1) to help decipher the difference between the two assessed accessions during the grain filling process. The unique and common DEGs between YG and ZG were subsequently identified using an UpSetR plot. T1 had the highest number of DEGs, with 587 upregulated and 362 downregulated genes in YG compared with ZG. T3 had the lowest number of DEGs, with only 133 upregulated and 80 downregulated genes in YG compared with ZG (Figs. 4A and 4B). The upregulated DEGs were annotated to 50 GO terms: with 22, 15 and 13 of these falling into the BP, CC and MF subcategories, respectively. Cellular process, single-organism process, and metabolic process were the top three classes in the BP category. Cell, cell part, and organelle were the top three classes in the CC category. Binding, catalytic activity, and nucleic acid binding transcription factor activity were the top three classes in the MF category. Similarly, the downregulated DEGs were annotated to 55 GO terms: with 24, 18 and 13 of these falling into the BP, CC and MF subcategories, respectively (Fig. S3). Cellular process, metabolic process, and single-organism process were the top three classes in the BP category. Cell, cell part, and organelle were the top three classes in the CC category. Binding, catalytic activity, and transporter activity were the top three classes in the MF category. Figures 4C and 4D show the significantly enriched pathways of the DEGs in YG and ZG at each stage of the grain filling process.

DEGs related to ATP-binding cassette (ABC) transporters

The ‘ABC transporter’ pathway was significantly enriched in the upregulated DEGs in the YG accession at late stage of grain filling (T5 and T6). Seven DEGs were related to ABC transporters which in the YG accession had a significantly higher expression than they did in the ZG accession (Fig. 5). In YG, Seita.5G283000 (ABCB1-like) and Seita.7G040200 (ABCC1-like) showed a higher expression at the T1 and T5 stages of grain development; Seita.5G283500 (ABCB1-like) was significantly upregulated at the T5 and T6 stages of grain development; Seita.4G239600 (ABCC8-like) and Seita.7G039700 (ABCC1-like) had higher expression levels at the T5 stage of development; and Seita.5G189700 (ABCB1-like) and Seita.1G055700 (ABCB1-like) showed their highest degree of expression at the T6 stage of development.

DEGs related to the plant hormone signal transduction pathway

Upregulated DEGs in YG were enriched in the “plant hormone signal transduction” pathway at grain filling stages T2 and T4 (Fig. 6). Five DEGs related to abscisic acid (ABA) signal transduction were identified, including two PYR/PYL, one PP2C, one SnRK, and one ABF gene (Fig. 6A). These genes were highly expressed at most stages of the grain filling process of the YG accession, especially at T2 stage. Four DEGs related to the auxin (IAA) signal transduction pathway were identified, including one ARF, two SAUR and one GH3. Most of these genes were highly expressed in YG at different grain filling stages (Fig. 6B). Five jasmonic acid (JA) signal transduction-related DEGs (one JAR1 and four JAZ genes) were also identified. Compared with the ZG accession, JAR1 had a significantly higher expression during the early stages, however, the expression levels of 4 JAZ genes were lower in YG during the early stages but increased rapidly during the late stages (Fig. 6C). In addition, 4 DEGs related to salicylic acid (SA) signal transduction were upregulated at most grain filling stages in YG (Fig. 6D). Notably, TGA was strongly upregulated from grain filling stages T2 to T6 in YG compared to ZG. Moreover, brassinoid steroids (BR) and cytokinin (CK) signal transduction downstream genes had higher levels of expression at most stages in YG when compared to ZG (Figs. 6E, 6F).

DEGs related to the GA biosynthesis pathway

The Diterpenoid synthesis pathway is a key pathway for GA biosynthesis. This pathway was significantly enriched in the upregulated DEGs in YG at both the preparation and initial stage of grain filling (T1 and T2). Seven gene families consisting of 12 genes involved in GA biosynthesis were identified as part of this analysis (Fig. 6G). GA biosynthesis genes in YG had a significantly higher expression than they did in ZG. In YG, the CPS, KS, KAO, GA20OX, GA2OX, and GA3OX encoded by loci Seita.9G351400, Seita.1G203500, Seita.6G089500, Seita.2G045400, Seita.5G147400, and Seita.3G103700, respectively, were all significantly upregulated at the T1 stage of grain filling. The CPS transcript encoded by loci Seita.1G203600 was significantly upregulated at grain filling stages T1 and T2. Similarly, two KO, one CPS, and one KS gene, encoded by Seita.4G235500, Seita.4G235600, Seita.2G144900, and Seita.2G144400, respectively, were significantly upregulated at the T1 and T4 stages of grain filling. The GA2OX encoded by Seita.3G148100 loci were significantly upregulated at grain filling stages T2 and T6.

DEGs related to photosynthesis and the photosynthesis-antenna protein pathway

The expression levels of genes involved in photosynthesis decreased during grain filling in both accessions (Fig. 2C). Besides, DEGs associated with photosynthesis and photosynthesis-antenna protein which in the YG accession had a significantly lower expression than they did in the ZG accession at the middle and late stages (Fig. 7).

We identified 29 DEGs associated with photosynthesis: nine involved in photosystem I, 12 in photosystem II, one in the cytochrome b6f complex, four in photosynthetic electron transport, and three associated with F-type ATPase function (Fig. 7). Thirteen DEGs associated with the photosynthesis-antenna protein were identified: seven related to LHC I and six to LHC II. Notably, all the DEGs associated with photosynthesis and the photosynthesis-antenna protein were highly expressed in the ZG accession, compared to the degree of their expression in the YG accession, especially during the middle and late stages. This phenomenon suggested that ZG could maintain a relatively high kernel photosynthetic capacity during grain filling compared to YG.

qRT-PCR verification of RNA-seq gene expression data

A total of 13 randomly selected DEGs (Seita.3G020300, Seita.7G193900, Seita.1G251000, Seita.7G001200, Seita.4G210500, Seita.8G118500, Seita.8G055600, Seita.9G022400, Seita.7G291500, Seita.8G051500, Seita.4G137900, Seita.9G346800 and Seita.8G199700) were used to confirm the validity of the sequencing data. The results indicated that the expression pattern of all genes in the qRT-PCR assays were essentially consistent with the sequencing data, indicating that the RNA-seq analysis was reliable (Fig. 8).