Two peanut cultivars exhibiting different stress tolerance under low calcium conditions
To investigate the phenotypic responses of peanuts under different calcium fertilization, plant morphological and agronomic traits were evaluated. As shown in Fig. 1A-1H, there were significant differences in phenotypic changes between “Lanshan” and “XH2008” when calcium shortage occurred. In terms of pod fulfilling, “XH2008” showed high rate of empty and bad pods while the phenotypic change in economic traits was not obvious in “Lanshan”. This can also be observed in Fig. 1I-1M, where full pod rate, full pod weight and yield per plant didn’t show significant changes in “Lanshan” while “XH2008” showed significant decrease under low-calcium condition (Fig. 1I-1K). As for plant architecture traits such as main stem height and lateral branch length, both cultivars were significantly decreased after the application of calcium fertilization (Fig. 1L-1M). Among all the phenotypic traits collected, only branch number neither showed significant genotypic effects nor genotypic-by-treatment interactions (Table 1). The correlation analysis among all the phenotypic traits was shown in Fig. 1N.
Table 1
The two-way ANOVA for seven phenotypic traits with the significance of genotypic and treatment main effects and their interactions.
Phenotype
|
Genotype
|
Treatment
|
Genotype × Treatment
|
Branch Number
|
0.500
|
***
|
0.067
|
Lateral Branch Length
|
****
|
****
|
*
|
Main Stem Height
|
****
|
****
|
***
|
Full Pod Rate
|
****
|
****
|
****
|
Full Pod Weight
|
****
|
****
|
****
|
Pod Setting Rate
|
**
|
*
|
*
|
Yield per Plant
|
***
|
*
|
*
|
Note: Genotypes included “XH2008” and “Lanshan”; Treatment levels included calcium deficiency and sufficiency levels; “*” represents P < 0.05, “**” represents P < 0.01, “***” represents P < 0.001, “****” represents P < 0.0001. |
Transcriptome profiles from two peanut cultivars treated under different calcium nutrition
In order to investigate differences in molecular mechanisms that lead to different responses to low calcium in two cultivars, the treated leaves of “Lanshan” and “XH2008” were sequenced. To briefly summarize, a total of approximately 547.94 million raw reads were generated from 12 cDNA libraries by RNA-Seq. After deleting adapter sequences as well as filtering low-quality and N-containing reads, 100% raw reads were confirmed as clean reads (Table S1). A total of 547 million clean reads were filtered. The percentage of clean reads mapped to the reference genome arahy.Tifrunner.gnm1.KYV3 ranged from 75.93 to 82.75% (Table S1), suggesting that the transcriptome sequencing quality was sufficient for further analyses.
Differentially expressed genes (DEG) analysis, enrichment analysis and qRT-PCR validation
Generally, a threshold of |log2FC| ≥ 2 and FDR < 0.05 was used to screen out differentially expressed genes (DEGs). The number of DEGs after filtration were 485, 58 and 61, respectively (NoCal vs Cal, S-NoCal vs S-Cal, and L-NoCal vs L-Cal) (Fig. 2A). Among these DEGs, 436 up-regulated and 168 down-regulated DEGs were identified. Of these DEGs between NoCal and Cal, 377 were up-regulated and 108 were down-regulated (Fig. 2A). Figure 2A identified 47 up-regulated and 11 down-regulated DEGs in small-seed cultivar while 12 up-regulated and 49 down-regulated DEGs in large-seed cultivar under calcium sufficiency when compared to calcium deficiency condition. On the whole, the number of up-regulated DEGs is higher than down-regulated DEGs, except large-seed cultivar “XH2008”. The Venn plots suggested that only one gene (LOC112797634) was common DEGs among the three comparison groups (Fig. 2B). Only 3 DEGs shared between “S” and “L” treatment groups (Fig. 2C), and the gene expression heatmap was shown in Fig. 2D. Among these 3 DEGs, LOC112766976 was uncharacterized with no function annotation. LOC112776736 and LOC112797634 were annotated as phosphoenolpyruvate carboxylase kinase 1 and UDP-glycosyltransferase 74G1, respectively. Phosphoenolpyruvate carboxylase kinase 1 belongs to a calcium-independent kinase, which is involved in light-dependent phosphoenolpyruvate carboxylase phosphorylation. UDP-glycosyltransferase 74G1 is involved in steviol glycosides biosynthesis.
To experimentally confirm the results of RNA-Seq data, six genes were selected to perform qRT-PCR (Table S2). As shown in Fig. 3, the selected genes had consistent expression patterns between RNA-Seq and qRT-PCR. This signifies that the RNA-seq data was of high-quality.
Co-expression network analysis identified key modules correlated with low-calcium tolerance
To identify the expression of genes related to low-calcium stress in two peanut cultivars, a gene co-expression network was constructed using WGCNA. The 2,650 selected genes were assigned to ten merged co-expression modules (Figure S1). As shown in Figure S2, we successfully identified two modules significantly associated with two peanut cultivars under different calcium treatments. The green module was only positively correlated with “Lanshan” under calcium fertilization, while the blue module was only positively correlated with “XH2008” without calcium fertilization. These two modules were possibly related to different responses of two cultivars under different calcium treatments. Even though brown and yellow modules showed differences between two cultivars, the changes were not found between different calcium treatments, and therefore, these two modules were not identified as interested ones throughout the current study. Additionally, a module-trait relationships analysis was performed using module eigengenes and phenotypic data (Table S3). As shown in Fig. 4, the MEgreen module was positively correlated with full pod rate (r = 0.54, p = 0.03) and full pod weight (r = 0.52, p = 0.04). In contrast, the MEblue module had a significant negative correlation with yield per plant (r = -0.75, p = 8e-04), full pod rate (r=-0.79, p = 3e-04), and full pod weight (r=-0.77, p = 5e-04), while positively correlated with branch number (r = 0.57, p = 0.02). The identification of peanut genotype-specific modules under calcium deficiency was particularly important. Based on the above results, the MEgreen and MEblue module were related to different responses in small- and large-seed cultivars respectively under calcium deficiency.
GO and KEGG enrichment analysis of the key modules
Within the MEblue module, we screened 514 genes to perform KEGG analysis (Fig. 5A). Multiple crucial pathways involved specifically in “XH2008” were determined, which included “plant-pathogen interaction”, “flavonoid biosynthesis”,
“isoflavonoid biosynthesis”, “phenylpropanoid biosynthesis”, “circadian rhythm”, “MAPK signaling pathway” and “ABC transporters”. All of the 61 DEGs between L-NoCal and L-Cal groups were covered by MEblue module and were further subjected for KEGG enrichment analysis. The major KEGG pathways enriched by these 61 DEGs were “glycosphingolipid biosynthesis” and “flavonoid biosynthesis” (Fig. 5B). The relative KEGG pathways of DEGs were summarized in Table S4. Within the MEgreen module, we screened 239 genes to perform KEGG analysis (Fig. 5C). Multiple crucial pathways such as “glycerophospholipid metabolism”, “glycerolipid metabolism”, “vitamin B6 metabolism”, “riboflavin metabolism”, “glycolysis/gluconeogenesis” and “ether lipid metabolism” were involved specifically in “Lanshan” cultivar under different calcium treatment. All of the 58 DEGs between S-NoCal and S-Cal groups were covered by MEgreen module and were further subjected in KEGG enrichment analysis. The major KEGG pathways enriched by these 58 DEGs were involved in “alpha-linolenic acid metabolism”, “tryptophan metabolism”, and “glyoxylate and dicarboxylate metabolism” (Fig. 5D). The relative KEGG pathways of DEGs were summarized in Table S5.
GO enrichment analyses including three categories were also conducted in these DEGs. Figure 6A-6C and Fig. 6D-6F showed GO term enrichment results of DEGs specifically in “XH2008” and “Lanshan”, respectively. In terms of cellular component, “XH2008” was mainly involved in mitochondria and photosystem, while “Lanshan” mainly involved in secondary plasmodesma, sieve plate and symplast. In terms of biological process, “XH2008” was mainly enriched in spermidine hydroxycinnamate conjugate biosynthetic process, response to chitin, fruit development and seed maturation, while “Lanshan” was mainly enriched in response to cycloheximide, cellular response to phosphate starvation, peptidyl-serine phosphorylation and intracellular signal transduction.
Hub genes involved in calcium deficiency tolerance via WGCNA
The key genes associated with low-calcium tolerance in blue and green modules were selected using WGCNA. As a result, a total of 128 and 35 genes were selected as hub-gene sets in blue and green modules, respectively. Gene network analysis was further conducted by Cytoscape software. Based on the degree value of each node, top 20 genes with connectivity degree larger than 97 and 17 were regarded as “hub genes” in blue and green modules, respectively (Fig. 7A-7B). KEGG enrichment analyses were separately conducted for the above two hub-gene sets (Fig. 7C-7D). Gene annotation and network analysis information of the hub-gene sets were provided in Table S6 & S7.
Among these hub-gene sets, 7 and 1 hub genes were identified as differentially expressed hub genes within blue and green module, respectively (Fig. 7E & Table 2). Among the seven hub-DEGs identified in blue module, LOC112755784 and LOC112755200 were characterized as probably WRKY transcription factors 40 and 33, respectively, which encode WRKY DNA-binding domains and play key roles in plant-pathogen interaction as well as mitogen-activated protein kinase (MAPK) signaling pathway. MAPK signaling pathway is known as sensing the extracellular stimuli and conduct signaling transduction process. LOC112766362 and LOC112764075 were characterized as NAD(P)H-dependent 6'-deoxychalcone synthase, belonging to keto reductase family, and related to flavonoid biosynthesis pathway. LOC112743667 encodes malonyl-CoA:anthocyanidin 5-O-glucoside-6''-O-malonyltransferase, which is related to isoflavonoid biosynthesis pathway. LOC112784753 was annotated as UDP-glycosyltransferase 83A1-like, belonging to UDP-glucoronosyl and UDP-glucosyl transferase family. The differential expression pattern of these 7 hub genes was shown in Fig. 7F. The molecular mechanisms related to the generation of secondary metabolites in the large-seed cultivar was proposed (Fig. 8) and the gene expression modes in the relevant pathways were provided in Table S8. The only one hub-DEG in green module was LOC112776736, which encodes phosphoenolpyruvate carboxylase kinase 1, which is a Ca2+-independent enzyme activated by a process involving protein synthesis in response to a range of signals in different plant tissues.
Table 2
Gene annotation and KEGG pathway for hub-DEGs in blue and green modules
Gene ID
|
Module
|
GeneBank Annotation
|
KEGG Pathway
|
LOC112784753
|
blue
|
UDP-glycosyltransferase 83A1-like
|
NA
|
LOC112766362
|
blue
|
NAD(P)H-dependent 6'-deoxychalcone synthase
|
Flavonoid biosynthesis
|
LOC112764075
|
blue
|
NAD(P)H-dependent 6'-deoxychalcone synthase
|
NA
|
LOC112755784
|
blue
|
Probable WRKY transcription factor 40
|
Plant-pathogen interaction
|
LOC112755200
|
blue
|
Probable WRKY transcription factor 33
|
MAPK signaling pathway;
Plant-pathogen interaction
|
LOC112743667
|
blue
|
Malonyl-CoA:anthocyanidin 5-O-glucoside-6''-O-malonyltransferase
|
Isoflavonoid biosynthesis;
Flavone and flavonol biosynthesis
|
LOC112705997
|
blue
|
Uncharacterized
|
NA
|
LOC112776736
|
green
|
Phosphoenolpyruvate carboxylase kinase 1
|
NA
|