Comparative transcriptome analysis during seeds development between two soybean cultivars

Plant materials

Soybean plants of the cultivars A7 (Yudou 12) and A35 (Fendou 53) were grown and total RNA was extracted from seeds as described in our previous research (Zhao et al., 2018b). For each stage, three independent RNA extractions were performed and used in RNA-seq and qRT-PCR experiments.

Determination of protein and oil contents of seeds

Determination of the oil content in soybean: the mature seeds were ground to powder which was transferred into 10 mL glass tubes. Oil was extracted with ligarine and total lipids (TL) were determined (Dong et al., 2011). The oil content of the soybeans was calculated on the basis of dry weight of the seeds. The content of proteins was determined by Kjeldahl’s method. One gram of fine seed powder was used for digestion, cooling and distillation-titration. The percentage of bean protein was the total nitrogen percentage multiplied by 6.25. Each experiment was performed three times.

cDNA library construction and transcriptome sequencing

The total RNA (2 μg) was sent to GENEWIZ (SuZhou, China) for sequencing, assembly and clustering analysis. The total RNA was quantified and qualified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and a NanoDrop system (Thermo Fisher Scientific Inc., Waltham, MA, USA). One μg total RNA with RIN value above 7 was used for cDNA library preparation. First, double-strand cDNA was synthesized, and then treated to repair both ends and add a dA-tailing, followed by a T-A ligation to add adaptors to both ends. Size selection of adaptor-ligated DNA was then performed using AxyPrep Mag PCR Clean-up (Axygen, Union City, CA, USA), and fragments of ~360 bp (with the approximate insert size of 300 bp) were recovered. Each sample was then amplified by PCR for 11 cycles using P5 and P7 primers, with both primers carrying sequences which can anneal to the flow cell to perform bridge PCR, and with the P7 primer carrying a six-base index allowing for multiplexing. The PCR products were cleaned up using AxyPrep Mag PCR Clean-up (Axygen, Union City, CA, USA), validated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), and quantified by Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA).

Then libraries with different indices were multiplexed and loaded on an Illumina HiSeq instrument according to manufacturer’s instructions (Illumina, San Diego, CA, USA). Sequencing was carried out using a 2 × 150 bp paired-end (PE) configuration. The raw data has been deposited in NCBI with the following accession numbers: A35-2: SRR12283471, SRR13090079, SRR13090078; A35-4: SRR12283470, SRR13090077, SRR13090076; A35-6: SRR12283469; A7-2: SRR12283468, SRR13090082, SRR13090083; A7-4: SRR12283467, SRR13090081, SRR13090080; A7-6: SRR12283466, SRR13100636, SRR13100635.

RNA-seq data analysis

The software of Trimmomatic v0.30 (http://usadellab.org/cms/index.php?page=trimmomatic) was used to remove the technical sequences for the high quality clean data (Bolger, Lohse & Usadel, 2014). As for mapping, the reference genome sequences were downloaded from Phytozome (https://phytozome-next.jgi.doe.gov/info/Gmax_Wm82_a2_v1, version Glyma 2.0, 975 Mb), and then the clean data were aligned to the reference genome using Hisat2 (v2.0.1) (http://daehwankimlab.github.io/hisat2/). Differential expression analysis used the DESeq Bioconductor package, a model based on the negative binomial distribution (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html) (Anders & Huber, 2010, 2012). After controlling for the false discovery rate (Benjamini & Hochberg, 1995), the P-value threshold of <0.05 was used for detecting significant differential gene expression. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed according to the methods described in Xin et al. (2019).

Principal component analysis

Principal component analysis was performed to visualize sample-to-sample distances, as described by Lee et al. (2014). This analysis was conducted in R (R Core Team, 2017).

Quantitative real time-PCR analysis and statistical tests

The methods for Quantitative Real Time-PCR (qRT-PCR) analysis and statistical tests were performed as described previously (Zhao et al., 2018b). The specific primers are shown in Table S1.

Phenotype identification and transcriptome sequencing

A7 and A35 are two cultivars of soybean from China. Their oil and protein contents in mature seeds vary greatly, with the average protein content of A7 (50.8%) being significantly higher than in A35 (38.1%), while the oil content in A7 (17.6%) is much lower than in A35 (22.11%) (Fig. 1B). To study the global changes in gene expression during seed development, seeds from both cultivars were collected and analyzed by transcriptome sequencing at different timepoints (weeks after flowering or WAF) resulting in the following samples: A7-2 WAF, A7-4 WAF, A7-6 WAF, A35-2 WAF, A35-4 WAF, A35-6 WAF (Fig. 1A).

All the sequencing data are shown in Table 1. A total of 139 and 153 million clean reads were acquired in A7 and A35, respectively. The reads were highly matched to the soybean reference genome, and the statistical results showed that 54,665,387, 44,537,617, 40,763,561, 39,120,215, 44,375,847 and 41,099,599 reads were mapped for A7-2, A7-4, A7-6, A35-2, A35-4 and A35-6, respectively, with 90.3% average matching rate. The saturation and evenness of the transcriptomes were further analyzed to estimate the sequencing depth and evenness. According to the gene expression levels, the saturation curves were analyzed in four classifications: below 25%, 25–50%, 50–75% and 75–100% (Fig. S1). The expression evenness was analyzed through the whole transcripts from 5′ to 3′ (Fig. S2). The results showed that the overall quality of sequencing in this study was high and covered the vast majority of expressed genes. In addition, the correlation analysis of different samples showed that the average correlation of all the comparison pairs was more than 90% (Fig. S3).

DEGs in paired comparisons of A7 and A35

In order to detect significant DEGs, the P-value was adjusted as above, and the log2 foldchange > 2 and P-value < 0.05 were set as thresholds (Fig. S4). In total, the DEGs from different developmental stages and different soybean cultivars have been systematically compared and valued. Comparing different time points (WAF) in the cultivar A7, 12,638 DEGs were found in A7-2-VS-A7-4 (7,237 up-regulated and 5,401 down-regulated); 7,595 DEGs were found in A7-4-VS-A7-6 (4,196 up-regulated and 3,399 down-regulated); and 17,035 DEGs were found in A7-2-VS-A7-6 (9,736 up-regulated and 7,299 down-regulated) (Fig. 2A). Similarly, in A35 the number of DEGs in A35-2-VS-A35-4 was 9,364 (6,134 up-regulated and 3,230 down-regulated); in A35-4-VS-A35-6 it was 7,375 DEGs (4,034 up-regulated and 3,341 down-regulated); in A35-2-VS-A35-6, it was 14,231 DEGs (8,828 up-regulated and 5,403 down-regulated). Generally, the DEGs in A7 and A35 mainly occurred in the 2- and 4- WAF pair and in the 2- and 6-WAF pair, and many more up-regulated genes were observed than down-regulated genes, especially in the early developmental stages (Fig. 2). In seeds, the storage nutrients gradually accumulated from the early stage of the seed development (two WAF) to the late stage of seed development (six WAF). On this basis, we conclude that a large number of genes are activated mainly around two WAF to meet the needs of energy and materials during the accumulation of storage nutrients in seeds.

Furthermore, the comparisons between A7 and A35 at the same stages were also performed. Overall, 1312 DEGs (725 up and 587 down) were found in both A7 and A35 (Fig. 2B). However, at different developmental stages a divergence of the DEGs was observed. Firstly, comparing A7 and A35, the DEGs mainly occurred in 2-WAF seeds, with the number of DEGs in 2-WAF seeds (4,818) being much higher than in the 4- (1,371) and 6-WAF seeds (1,515). This observation was supported by the results of the principal component analysis (PCA) showing a major difference between A7 and A35 at the early stage of seed development (two-WAF seeds) (Fig. S5). Secondly, in the 2-WAF seeds, the number of upregulated DEGs (3,317) was significantly higher than the number of downregulated genes (1,501) in A35 compared to A7, which indicated more genes were activated in A35 than A7 during the early developmental stages of seeds. In addition, the time point specific DEGs were also compared, and the results were consistent with above ones, in which 3,806, 446 and 708 specific DEGs have been identified at two WAF, four WAF and six WAF, respectively.

Analysis of DEG clustering and pathway enrichment

To better understand the trends of DEG regulation during the seed development, hierarchical clustering was performed and presented in a boxplot dendrogram (Fig. S6). The dendrogram provided a clear overview of the clade structure, in which the DEGs clustered to the up-regulated or down-regulated clades, according to their expression trends. In the heatmap (Fig. S6), 25,139 DEGs were shown, but only half of them (13,262) were annotated in GO involving biological process (5,631), cellular component (1,200) and molecular function (6,431). Furthermore, the DEGs concentrated in the functional groups such as catalytic activity, binding, metabolic process, cellular process and biological regulation, which were consistent with the processes of seed development including energy metabolism, nutrients accumulation and cell proliferation (Fig. 3A; Fig. S7).

To further find out the involved metabolic pathways of the DEGs, the data were mapped against the KEGG pathway database (Table S5). In A7, a total of 85 pathways with Q value < 0.05 were identified. Fifty eight, 46 and 33 pathways were mapped in the A7-2-VS-A7-4, A7-2-VS-A7-6 and A7-4-VS-A7-6 paired comparisons, respectively. In A35, 78 pathways were totally identified: 46 pathways were in A35-2-VS-A35-4, 37 in A35-2-VS-A35-6 and 28 in A35-4-VS-A35-6, respectively. Notably, eight and four significantly mapped pathways were common to the three paired comparisons in A7 and A35, respectively (Fig. 3B). In A7, the common pathways included ko05203 (Viral carcinogenesis), ko01524 (Platinum drug resistance), ko00073 (Cutin, suberine and wax biosynthesis), ko05418 (Fluid shear stress and atherosclerosis), ko00480 (Glutathione metabolism), ko05034 (Alcoholism), ko04110 (Cell cycle) and ko04914 (Progesterone-mediated oocyte maturation), which were mainly involved in cell cycle and stresses responses. However, in A35, the common pathways were ko00052 (Galactose metabolism), ko00073 (Cutin, suberine and wax biosynthesis), ko00591 (Linoleic acid metabolism) and ko00500 (Starch and sucrose metabolism), mainly involved in sugar and fatty acids metabolism. Only one common pathway, namely ko00073, occurred in both A7 and A35, which was related to the biosynthesis of cutin, suberine and wax in cell wall (Fig. 3). When specific developmental stages of A7 were compared with A35, lipid-related pathways (ko00592, ko00062 and ko00071) and sugar-related pathways (ko00196, ko00195, ko00710 and ko00010) were found at the 2-WAF stage, while amino acid metabolism pathways (ko00360 and ko00270) were found at 4- and 6-WAF stages, respectively (Table S5). The divergence of the pathways in different developmental stages of soybeans showed conformity with the differences in nutrient contents between A7 and A35.

Exploration of DEGs that may contribute to the oil and SSP contents in soybeans

Given the significant difference in oil and protein contents in A7 and A35, the DEGs involved in the synthetic and metabolic pathways of lipid, fatty acids, amino acids and proteins during the seed development were specifically investigated (Fig. 4). Firstly, the change trends of up- and down-regulated DEGs involved in the lipid metabolism were V-shaped in A7 and A35 from 2-VS-4, 4-VS-6 to 2-VS-6. Correlation analysis showed that the trends were consistent between A7 and A35 (Fig. 4A). Similarly, the change trends of DEGs related to the amino acid metabolism were also V-shaped in A7 and A35 (Fig. 4D). On the other hand, the change trends of DEGs involved in fatty acid and protein metabolism were different in A7 and A35. As for the fatty acids, the numbers of up- and down-regulated DEGs in A7 were similar in 2-VS-4 and 4-VS-6 pairs, but sharply decreased in 2-VS-6. However, in A35 the number of down-regulated DEGs was increased linearly from 2-VS-4, 4-VS-6 to 2-VS-6, while the trend of the up-regulated DEGs was V-shaped (Fig. 4B). Notably, the trends of DEGs involved in protein synthesis were opposite in A7 and A35. In A7, the number of up- and down-regulated DEGs was first increased and then decreased, while in A35 it was the other way around (Fig. 4C). Notably, the divergent trends of numbers of DEGs involved in fatty acid and protein synthesis were consistent with the difference in oil and protein contents in A7 and A35, which indicated that these DEGs probably contribute to the content differences in A7 and A35.

To further reveal DEGs involved in oil and protein content differences in A7 and A35, the significant DEGs were matched to the QTL database in SoyBase (https://www.soybase.org/). In this way, two DEGs, Glyma.06G214800 and Glyma.07G034900, were found to be involved in oil and linoleic acid contents. Glyma.06G214800 encodes 3-ketoacyl-CoA synthase (KCS) catalyzing the elongation of C18 fatty acids while Glyma.07G034900 encodes linoleate 9S-lipoxygenase (LOX1-5) which generates a peroxide by adding oxygen in the double bond of linoleic acid. In our study, the expression levels of Glyma.06G214800 and Glyma.07G034900 were down-regulated during the seed development in A35. However, the expression of Glyma.06G214800 was not detected in A7, while the expression levels of Glyma.07G034900 in A7 were significantly higher than in A35, which probably is related to the difference in oil content between A7 and A35 (Fig. 5A; Table S6).

No specific DEGs involved in protein content were matched to the QTL database. However, previous studies have revealed that the expression of SSPs is mainly controlled by transcription factors during the seed filling stages, the master regulators include LEC1, Leafy Cotyledon 2 (LEC2), Abscisic acid-intensitive 3 (ABI3) and FUSCA3 (FUS3) (reviewed by Fedorova et al., 2002). In our study, the expression of LEC1s, LEC2s, FUS3s and ABI3s showed some divergence, in which the expression trends of LEC1s (Glyma.07G268100 and Glyma.17G005600) and ABI3 (Glyma.02G099500, Glyma.01G087500, and Glyma.10G204400) were basically consistent in A7 and A35, while LEC2s (Glyma.20G035800 and Glyma.20G035700) and ABI3 (Glyma.20G186200) expression was different (Fig. 5; Table S6). Most of the consistently expressed transcription factors (TF) were higher in A35 than A7 except for ABI3 (Glyma.02G099500 and Glyma.01G087500) (Fig. 5; Table S6). Considering the positive relationship of TF expression with the content of SSPs, the two DEGs, ABI3s (Glyma.02G099500 and Glyma.01G087500) and LEC2s (Glyma.20G035800 and Glyma.20G035700), might be correlated with the difference of SSPs between A7 and A35. Furthermore, a sugar transporter, SWEET10a (Glyma.15G049200), contributing to the increase in seed weight and oil content and reduced protein content was identified by the QTL database search. Its expression patterns were different in A7 and A35, with an inverted V shape in A7 and linear decline in A35. The biggest divergence between A7 and A35 occurred at the 2-WAF stage, where the expression in A35 is much higher than in A7 (p < 0.001, Fig. 5; Table S6).

Expression flux in the Fatty acids synthesis pathway

The pathways of fatty acid metabolism mainly include fatty acid biosynthesis (ko00061), fatty acid elongation (ko00062) and fatty acid degradation (ko00071) (Fig. 6; Fig. S8). All the DEGs and their expression trends in these pathways have been marked in Fig. 6 to identify their expression flux. The highest expression levels in the whole pathway were observed for fatty acid dehydrogenase 2 (FAD2s) encoded by Glyma.10G278000 and Glyma.20G111000. Both DEGs were specifically expressed in seeds and the expression levels were increased during the seed development. The expression levels of FAD2s were higher in A35-2WAF than in A7-2WAF, which suggested that FAD2 genes might be essential for the content of oil in A7 and A35. Notably, the expression levels of key enzymes of the fatty acid elongation pathway were drastically reduced during the development of seeds in A7 and A35 (Fig. 6) and their expression levels in A35 were lower than in A7. The results might explain why the fatty acids in soybean seeds were mainly present in the form of C18. As for fatty acid degradation, the expression of enoyl-CoA hydratase/3-hydroxyacyl-CoA dehydrogenase (MFP2/HAD) and acetyl-CoA C-acetyltransferase (atoB) was stable in the development of seeds in A7 and A35 (Fig. 6).

The qRT-PCR analysis of DEGs in soybeans

To confirm the RNA-seq results in our research, we selected 6 significant DEGs related to the fatty acid metabolism and the synthesis of amino acids for qRT-PCR analysis (Fig. 7; Table S6). In all, the expression results of qRT-PCR were consistent with the RNA-seq results. Among them, the expression of Glyma.13G035200 was significantly increased along the seed development, while Glyma.16G147300 was gradually decreased. The expression trend of Glyma.06G211300 was decreased from 2-WAF to 4-WAF seeds, and then was increased in 6-WAF seeds. The expression trends of Glyma.17G047000, Glyma.17G027600 and Glyma.06G183900, on the other hand, were bell-shaped, increasing from 2-WAF to 4-WAF seeds and then decreased in 6-WAF seeds. The results verified that our RNA-seq data were reliable.