Flowering-associated gene expression and metabolic characteristics in adzuki bean (Vigna angularis L.) with different short-day induction periods

Experimental material

Portions of this text were previously published as part of a preprint (https://doi.org/10.21203/rs.3.rs-3362672/v1). A late-maturity variety ‘Tangshan hong xiao dou,’ which is sensitive to short days, was selected as the experimental material. Plants were provided by the Adzuki Bean Breeding Research Group of the Institute of Grain and Oil Crops, Hebei Academy of Agricultural and Forestry Sciences (China). The experiment was conducted in the Teaching and Experimental Base of Hebei Agricultural University, Baoding (longitude: 115°47′, latitude: 38°87′) in 2022. As a typical short-day crop, adzuki bean is highly sensitive to short-day induction; thus, the flowering time and maturity stage are significantly earlier than those under a natural photoperiod.

Soil fertility

The experimental field was comprised of loam soil, and 400 g of compound fertilizer (N:P₂O₅:K₂O = 24:4:8) was applied to each plot (5-m long, 1-m wide) prior to sowing. The nutrient content of the experimental plot was measured after ploughing. The fertility of the cultivated soil layer (0–20 cm) of the experimental plot is shown in Table 1.

Environmental characteristics of the plot after shading

The environmental characteristics of the plot after shading treatment are shown in Table 2. Compared with the atmospheric environment, the light intensity of the plot after shading was near zero, the relative humidity was significantly increased by 21.41%, and the CO₂ concentration and temperature were slightly, but not significantly increased.

Short-day treatment

Compound fertilizer (400 g) was applied to the experimental plot, followed by ploughing, and then sowing which was carried out on June 24, 2022. Planting was performed in two rows, plant spacing and row spacing were 15 × 40 cm. Plants were grown under natural light until the true leaf unfolded, after which they were placed on a stainless-steel shelf (6-m long, 1.5-m wide, 1-m high) on the upper side of the plots and covered with an opaque cloth. The opaque cloth (10-m long, 5-m wide) comprised two layers of black cloth and two layers red cloth stitched together. The shading treatment involved 10 h of light and 14 h of dark with shading for 5, 10, or 15 d. By adopting a randomized block design, this experiment involved initiating shading treatment with the opaque cloth at 18:00 every day until 08:00 the following morning. After the shading treatment, all plants were grown to maturity under natural light. At the initial flowering and pod-setting stages, 0.4% potassium dihydrogen phosphate was sprayed and irrigation. There were three replicates for each shading treatment (SD-5d-1, SD-5d-2, and SD-5d-3; SD-10d-1, SD-10d-2, and SD-10d-3; SD-15d-1, SD-15d-2, and SD-15d-3), with a total of nine plots. After 5, 10, and 15 d of shading, the middle leaf samples of the top trifoliate leaves were uniformly removed at 09:00 from 5–6 plants in each group. After sampling, the leaves were immediately frozen in liquid nitrogen and stored at −80 °C. RNA extraction and sequencing were performed 1 week later.

Meteorological factor determination

Light intensity was measured 20–30 cm above the canopy of the adzuki bean community using a TES1332 illuminometer (provided by College of Plant Protection, Agricultural University of Hebei, TES1332; TES, Taipei, Taiwan). The CO₂ concentration was measured using an Li-6400 portable photosynthetic meter (LI-COR, Lincoln, NE, USA). Temperature and humidity were measured using a HOBO Pro V2 series data logger (U23-002; Onset Computer Corporation, Bourne, MA, USA), which automatically counted and recorded data every hour.

Growth index determination

Plant height, stem diameter, and leaf area were measured in three representative plants selected from each plot at the flowering, pod-setting, and grain-filling stages. Plant height was calculated as the distance from the true leaf to the growing point of the plant. The plant stem diameter at the true leaf was measured using a vernier caliper and calculated as: circumference = 2 * 3.14 * (diameter/2). Leaf area was measured using a YMJ-B leaf area measuring instrument (Hangzhou Huier Instrument Equipment Co., Ltd., Hangzhou, China).

Flowering characteristics

At the flowering stage, three representative plants were selected from each treatment. The advanced flowering days were recorded and the flowering promotion rate was calculated using the following formulae:

Advanced flowering days = number of days from emergence to flowering of plants in the control group − number of days from emergence to flowering of plants in the treatment group.

Flowering promoting rate (%) = ([number of days from emergence to flowering in the control group − number of days from emergence to flowering in the treatment group] × 100)/number of days from emergence to flowering in the control group.

Apical flower buds were removed after short-day induction. After sampling, the buds were fixed, dehydrated and made transparent, paraffin-embedded, and sectioned. The paraffin-embedded sections were stained with hematoxylin and eosin and then observed under a microscope (BS500/BS500-TR biological microscope). Flower bud differentiation in adzuki bean was divided into the pre-differentiation (PD), inflorescence primordium differentiation (IP), flower primordium differentiation (FP), sepal primordium differentiation (SP), petal primordium differentiation (PP), stamen and carpel primordium differentiation (SCP), and stamen and carpel structural differentiation (SCS) stages according to the division method described by Jin, Chen & Yu (1995).

RNA extraction and cDNA library construction

Total RNA was extracted using a mirVana miRNA Isolation Kit (Ambion, Austin, TX, USA) following the manufacturer’s protocol and DNA was digested using DNase. Eukaryotic mRNA was enriched using oligo coupled to magnetic beads, and the mRNA was fragmented by addition of interrupting reagents. The fragmented mRNA was used as the template for first-strand cDNA synthesis using random hexamers. A second-strand synthesis reaction system was used to synthesize second-strand cDNA, which was further purified using the kit. End repair was then performed, followed by A-tailing and sequencing adapter ligation, fragment size selection, polymerase chain reaction (PCR) amplification, and library construction. After qualification using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA), the constructed libraries were sequenced using a HiSeq2500 or HiSeq X Ten sequencer (Illumina, San Diego, CA, USA) to produce 125 or 150 bp paired-end reads, respectively.

Screening of DEGs

The fragments per kilobase of transcript per million mapped reads (FPKM) metric (Roberts et al., 2011), bowtie2 (Langmead & Salzberg, 2012), and cufflinks software (Roberts & Pachter, 2013) were used to analyze the transcript levels. Using the Vigna angularis reference genome (https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/190/045/) (Kang et al., 2015), the number of transcript (protein-coding) reads was obtained for all samples using cufflinks software. Genes with an average number of reads >2 were screened, and related data were standardized using the estimateSizeFactors function of the R package DESeq (Livak & Schmittgen, 2001). P-values and fold-change values of difference comparisons were calculated using the nbinomTest function. DEGs were identified based on a fold change ≥2 and a false discovery rate <0.05.

Functional annotation analysis of DEGs

Gene ontology (GO) annotations and functional classification of DEGs (P < 0.05 and fold change >2) were performed using Blast2 GO and WEGO. BLAST was used to align gene sequences to those in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database for annotation of biochemical pathways and identification of the regulatory‒metabolic network.

Quantitative reverse transcription–PCR validation

Under the same treatment, other adzuki bean samples were used for quantitative reverse transcription‒PCR (qRT‒PCR) analysis of eight selected DEGs to verify the accuracy of the sequencing results. RNA was extracted using an RNA extraction kit followed by reverse transcription to produce cDNA (HiScript II Q RT SuperMix for qPCR). Primer 5.0 was used to design the primers; the primer sequences are shown in Table 3. A QuantiFast^® SYBR^® Green PCR Kit was used to prepare the qRT‒PCR reaction system and the transcripts were detected on a fluorescence quantitative PCR instrument (CFX96; Harbin Deyuan Science and Technology Development Co., LTD, Harbin, China). The qRT‒PCR cycling conditions were as follows: predenaturation at 95 °C for 10 min, followed by 40 cycles of denaturation at 95 °C for 10 s, and annealing/extension at 60 °C for 30 s. Gene expression levels were calculated using the 2^−∆∆Ct method (Lin & Pang, 2019).

Data analysis

Raw data were filtered using NGS QC Toolkit software to obtain high-quality, clean reads. FPKM values were calculated using cufflinks software. All data were analyzed using analysis of variance, with three replicates, using Excel 2023 and SPSS 17.0 (SPSS Inc., Chicago, IL, USA). Duncan’s new multiple range test (at the 5% probability level) was used to test the differences among the mean values. Images were processed using BS500/BS500-TR biomicroscopy and Photoshop CS6.

Growth and flowering characteristics

Different short-day inducement periods had very different effects on apical flower bud differentiation and morphology, the number of days until early flowering, and the flowering promotion rate in adzuki bean (Fig. 1). Plant morphology differed under different short-day induction periods, with a longer short-day induction period resulting in a greater inhibitory effect on plant height of adzuki bean (Fig. 1B). We further observed apical flower bud differentiation and found that the apical flower bud was in the PD stage under the SD-5d treatment; whereas, it was in the SP stage under the SD-10d and SD-15d treatments (Fig. 1A). Plant height, stem diameter, and leaf area decreased with prolongation of short-day induction, and plant height showed significant differences among the three treatments at the flowering, podding and seed-filling stages (Fig. 2A). The stem diameter significantly decreased by 10.91% and 13.27% under the SD-15d treatment compared to that under the SD-5d treatment at the flowering and seed-filling stages, respectively; whereas, there was no significant difference in stem diameter between the SD-5d and SD-10d treatments. There were significant differences in stem diameter among the three treatments during the podding stage (Fig. 2B). Significant differences in leaf area were observed among the three treatments at all growth periods (Fig. 2C).

The flowering date was 12 and 7 d earlier in the SD-15d group than in the SD-5d and SD-10d groups, respectively (Fig. 1C), and the SD-15d group improved flowering promotion rates by 22.03% and 11.83% compared to the SD-5d and SD-10d groups, respectively (Fig. 1D). Advanced flowering occurred 6 d later in the SD-10d group than in the SD-5d group, and the flowering promotion rate was 10.20%. These results indicated that short-day induction promoted early flowering of adzuki bean by inhibiting vegetative growth. Short-day induction had a cumulative effect, as a longer short-day inducement resulted in an earlier flowering time and a higher flowering promotion rate.

Functional analysis of DEGs

Principal component analysis of the samples showed that the dispersion degree among individual samples was significant, and the first principal component (PC1) accounted for 62.71% of the variance. This separation underscored the reproducibility, accuracy, and overall reliability of the data (Fig. 3A). By comparing the data among the three groups using a volcano map, the total number of DEGs was 5,939, including 3,107 up-regulated genes and 2,832 down-regulated genes (Figs. 3B‒3D). When comparing SD-5d VS SD-10d, 2,044 DEGs were detected, including 961 up-regulated genes and 1,083 down-regulated genes (Fig. 3B). When comparing SD-5d VS SD-15d, 3,068 genes were detected, including 1,711 up-regulated genes and 1,357 down-regulated genes (Fig. 3C). When comparing SD-10d VS SD-15d, 827 genes were detected, including 435 up-regulated genes and 392 down-regulated genes (Fig. 3D). In addition, by comparing the three groups, 105 common genes were detected, of which 40 were up-regulated and 25 were down-regulated (Figs. 3E‒3G). In addition, the SD-5d VS SD-15d comparison showed the highest number of DEGs, indicating that a longer two short-day inducement treatment resulted in the detection of a higher number of DEGs.

Cluster analysis showed that the three treatments significantly affected gene expression (Fig. 3H). Among them, 40 genes were gradually up-regulated with extension of the short-day treatment time; whereas, 25 genes were gradually down-regulated. An additional 40 genes did not show obvious expression patterns, with 25 genes showing up-down-up regulation and 15 genes showing down-up-down regulation in all treatment groups. These results showed that 40 up-regulated and 25 down-regulated genes played key roles in adzuki bean growth and development; however, the other genes had no obvious roles, indicating that the function of different genes varied in response to different short-day induction times.

GO functional enrichment analysis

GO terms (DEGs > 2) were screened separately for biological processes, cell components, and molecular functions. GO functional enrichment analysis showed that light-related functions were found among DEGs from all three comparison groups. In the SD-5d VS SD-10d comparison, five up-regulated genes were related to photosystem I (PSI) and II (PSII). In addition, there was a circadian rhythm function related to light among the biological processes, which contained 14 genes, of which six were up-regulated and eight were down-regulated (Fig. 4A). The DEGs in the two comparison groups of SD-5d VS SD-15d and SD-10d VS SD-15d had light-driven functions in PSI and PSII; however, the number of genes involved differed. For the SD-5d VS SD-15d comparison, there were nine and 12 genes associated with PSI and PSII, respectively; the nine genes associated with PSI were all up-regulated, 11 of the 12 genes associated with PSII were up-regulated, and the remaining one gene was down-regulated (Fig. 4B). Both PSI and PSII functions contained five up-regulated genes in the SD-10d VS SD-15d comparison (Fig. 4C). These results suggested that the DEGS in the three comparison groups contained PSI and PSII functions, which are related to light (Fig. 4D), and that photosystem genes in leaf cells may have been activated after different short-day induction times. These DEGs were key regulators of adzuki bean seedling growth and development and were closely related to photoelectron transport under short-day-induced conditions.

KEGG pathway enrichment analysis of DEGs

KEGG pathway enrichment analysis of DEGs in the three groups was performed to select the top 20 signaling pathways with significant differences in each comparison group. Figure 5 shows that the DEGs identified in the SD-5d VS SD-10d and SD10d VS SD-15d comparisons were enriched in the antenna protein pathway and circadian rhythm pathway two light-related pathways, with different numbers of genes involved in these pathways. The antenna protein pathway in both comparison groups contained five genes that were all up-regulated. In the SD-5d VS SD-10d comparison, of the 16 DEGs in the circadian rhythm pathway, three were up-regulated and 13 were down-regulated (Fig. 5A). For the SD-10d VS SD-15d comparison, there were seven DEGs in the circadian rhythm pathway, including four up-regulated genes and three down-regulated genes (Fig. 5C).

Four light-related pathways were enriched in the SD-5d VS SD-15d comparison: antenna protein, circadian rhythm, plant hormone, and photosynthesis pathways, including 9, 17, 45, and 14 DEGs, respectively. The expression patterns of genes in different pathways differed. All nine DEGs in the antenna protein pathway were up-regulated. Of the 17 DEGs in the circadian rhythm pathway, three were up-regulated and 14 were down-regulated. A total of 45 DEGs were identified in the plant hormone signal transduction pathway, consisting of 27 up-regulated and 18 down-regulated genes. All 14 DEGs in the photosynthesis pathway were up-regulated (Fig. 5B). These results showed that all three comparison groups contained two pathways related to antenna proteins and circadian rhythms and that the SD-5d VS SD-10d and SD-5d VS SD-15d comparisons had the most significant degrees of enrichment in the circadian rhythm pathway (Fig. 5D). This suggested that different short-day induction periods, on the one hand, activated gene expression in the antenna protein pathway and promoted photosynthesis, and on the contrary, regulated photosynthesis, stomatal rhythmic opening and closing, leaf rhythmic opening and closing, and other important physiological processes, by regulating the biological clock. This may help adzuki bean plants to optimize their energy utilization efficiency and guarantee early flowering. Additionally, the circadian pathway played a key role in responding to different short-day induction periods.

Gene set enrichment analysis in the KEGG pathway

It is generally believed that a positive enrichment score (ES) value indicates that a certain functional gene set is enriched at the front of the ranked sequence and that the involved pathway is up-regulated, and a negative ES value indicates that a certain functional gene set is enriched at the rear of the ranked sequence and that the involved pathway is down-regulated (Subramanian et al., 2005). The core genes that appeared in the ES map indicated that this functional gene set had biological significance. Here, gene set enrichment analysis revealed that the antenna proteins and circadian rhythm pathways in the three comparison groups had the same enrichment scores (Fig. 6). The antenna protein pathway was up-regulated (Figs. 6A‒6C), the circadian rhythm pathway was down-regulated (Figs. 6D‒6F), and core genes were identified in both pathways. The core genes of the antenna protein pathway were enriched in the front of the sequencing sequence; whereas, those of the circadian rhythm pathway were enriched at the rear of the sequencing sequence in the three comparison groups, indicating that these core genes played important roles in the response of adzuki bean to short-day induction.

Heatmap cluster analysis was performed on the core genes of antenna protein and circadian rhythm metabolic pathways (Fig. 7). The circadian rhythm pathways in the SD-5d VS SD-10d, SD-5d VS SD-15d, and SD-10d VS SD-15d comparisons contained 12, 15, and 14 core genes, respectively, which were inversely down-regulated as the number of short-day induction days increased (Figs. 7A‒7C). The antenna protein pathways of the three comparisons, contained 13, 15, and 17 core genes, respectively, which were up-regulated as the number of short-day induction days increased (Figs. 7D‒7F). These results indicated that adzuki bean responded to different short-day induction periods by up-regulating key genes in the antenna protein pathway or down-regulating key genes in the circadian rhythm pathway, thus promoting early flowering.

Analysis of circadian and antenna protein pathway maps

The circadian rhythm involves a rhythmic change in plants in response to daylight length. Analysis of the circadian rhythm metabolic pathway showed that blue light promoted flowering and far-red light inhibited flowering in adzuki bean. After short-day induction, PHYB (i.e., a far-red light receptor interacted with PIF3) (i.e., a phytochrome interaction factor) to promote NDA transcription and production of LHY (i.e., a circadian clock regulatory protein). The PRR7 protein of the central oscillator-controlled internal clock inhibited PIF3 transcription and LHY production. In addition, PRR7 (i.e., a key gene in the photoperiodic flowering pathway) inhibited LHY protein expression, thus, inhibiting adzuki bean flowering. After short-day induction, one pathway regulating flowering involved the inhibition of complex formation between the morphogenesis regulator COP1 and the signal sensing protein SPA1 by the cryptochrome CRY-blue light receptor protein. This complex was down-regulated and ubiquitinated to inhibit transcription of HY5 (i.e., a negative regulator of photomorphogenesis), thus indirectly promoting photomorphogenesis. Another pathway regulating flowering was the CRY pathway. CRP is a blue light receptor protein that inhibited formation of a complex between COP1 and ELF3, thus inhibiting GI by ubiquitylation. As a negative regulator, GI down-regulated the function of the complex formed by FKF1 and GI; inhibited CDF1 binding to the CO promoter by ubiquitylation; and thus, initiated transcriptional activation of CO and up-regulation of FT. GI also directly promoted CO transcription, which can up-regulate FT, and indirectly promoted adzuki bean flowering (Fig. 8A). Figure 8B shows that changes in the expression levels of seven protein-coding genes (i.e., PIF3, PRR7, SPA1, HY5, ELF3, GI, and FKF1) differed among the SD-5d VS SD-10d, SD-5d VS SD-15d, and SD-10d VS SD-15d comparisons. More significant changes were observed in the SD-5d VS SD-15d comparison than in other comparisons, indicating that a longer short-day induction period caused greater gene expression change and that most genes were down-regulated.

Chloroplasts in green plants are comprised of PSI and PSII (Ferroni et al., 2016). PSI contains five light-trapping pigment protein complexes (i.e., Lhca1, Lhca2, Lhca3, Lhca4, and Lhca5) and PSII contains six (i.e., Lhcb1, Lhcb2, Lhcb3, Lhcb4, Lhcb5, and Lhcb6), of which Lhcb1–3 play major roles (Green, Pichersky & Kloppstech, 1991). Different short-day induction periods may affect the expression of genes encoding antenna proteins, allowing them to adapt to adzuki bean growth and development. Figures 9A and 9B show that the genes encoding the Lhca3, Lhcb1, Lhcb2, Lhcb3, Lhcb4, and Lhcb6 proteins were all up-regulated after a short-day induction period, and that there was a greater number of up-regulated genes encoding the Lhcb1, Lhcb2, and Lhcb3 proteins than of those encoding the other proteins. Figure 9B shows that the largest number of DEGs occurred in the SD-5d VS SD-15d comparison, indicating that the light signal reception ability was enhanced with an increase in the number of short-day induction days; that is, a longer short-day induction period resulted in stronger reception ability. These results indicated that adzuki bean responded to short-day induction by up-regulating PSI- and PSII-related genes. These genes may play key roles in regulating adzuki bean seedling growth and development under short-day induction conditions, and a longer induction period resulted in greater changes in gene expression levels.

Identification of DEGs by qRT–PCR

To verify the reliability and accuracy of the transcriptomic data, samples from the same treatment group were subjected to qRT‒PCR analysis. Eight genes with significant variability were screened from the circadian rhythm and antenna protein pathways, all of which were closely related to light-related functions and associated with flowering based on homologous sequence alignment (Table 4). qRT‒PCR was performed to verify the differential expression levels of these eight genes in SD-5d VS SD-10d, SD-5d VS SD-15d, and SD-10d VS SD-15d comparisons.

The LOC108331766 gene expression levels were significantly different in the SD-5d VS SD-10d and SD-5d VS SD-15d comparisons, but not in the SD-10d VS SD-15d comparison. The LOC108322606 gene showed significantly different expression levels in the SD-5d VS SD-10d and SD-10d VS SD-15d comparisons but not in the SD-5d VS SD-15d comparison. In the antenna protein pathway, LOC108345872 and LOC108344684 genes showed significantly different expression levels among the three groups. The LOC108335068 gene showed significantly different expression levels in the SD-5d VS SD-15d and SD-10d VS SD-15d comparisons, but not in the SD-5d VS SD-10d comparison. LOC108328079, LOC108333950, and LOC108338432 genes showed significantly different expression levels in the SD-5d VS SD-10d, SD-5d VS SD-15d, and SD-10d VS SD-15d comparisons, respectively, but not in the other comparisons. Therefore, the qRT‒PCR results for these eight genes were consistent with the RNA-seq results (Figs. 10A1B1, 10A2B2, 10A3B3, and 10A4B4). These eight genes in the circadian rhythm and antenna protein pathways showed significantly different expression levels, indicating that the transcriptome sequencing data were accurate and reliable.