The Arachis hypogaea “XiangHua 2008,” a peanut cultivar from the Hunan Agriculture University, was used as material. The red soil overburden including exchangeable calcium 0.74 cmol(1/2 Ca²⁺)/kg in Gengyuan practice base, Hunan Province of China, was used for the free Ca²⁺ deficient treatment. The soil fertilized with 100 kg plaster (CaO) per 667 m² was used as free Ca²⁺ sufficiency treatment. CaO was applied before sowing.

In peanut production, 15 days after pegging (DAP) was confirmed to be the most sensitive period for the free Ca²⁺ content in the soil. In this study, gynophores under free Ca²⁺-deficient treatment 1 (GD1), free Ca²⁺-deficient treatment 2 (GD2), free Ca²⁺-sufficient treatment 1 (GS1), free Ca²⁺-sufficient treatment 2 (GS2), as well as pods under free Ca²⁺-deficient treatment 1 (PD1), free Ca²⁺-deficient treatment 2 (PD2), free Ca²⁺-sufficient treatment 1 (PS1), and free Ca²⁺-sufficient treatment 2 (PS2), were harvested, frozen in liquid nitrogen and then stored at −70°C. The 8 samples were used to detect global changes in gene expression. Two biological replicates were used in this study. The mixed samples were used in transcriptome sequencing.

Total RNA was isolated from the harvested materials using the total RNA extraction Kit (TaKaRa, Inc., Dalian, China) according to the manufacturer's instructions. Agilent 2100 and NanoDrop were used to detect RNA quality and purity. Only high-quality RNA samples were chosen for RNA-seq analyses. Beads with Oligo (dT) were used to enrich for mRNA and the RNA was fragmented into short fragments. First-strand cDNA was synthesized using random hexamer primers and buffer, while dNTPs, RNaseH and DNA polymerase I were added to synthesize the second strand (Wang et al., 2009). After purification, the short fragments were then connected using sequencing adapters. Finally, the library was sequenced from directions on Illumina HiSeq™ 2000 System (Illumina, San Diego, CA). Data analysis and base calling were achieved by applying the Illumina instrument software. Whole dataset has been deposited in the NCBI Sequence Read Archive with accession number SRX1795063.

Differential expression profiling was conducted among the gynophores and pods under free- Ca²⁺-deficient or sufficient treatments. Prior to further analysis, raw data should be filtered to decrease data noise and obtain clean data (Cock et al., 2010). Gene expression level was quantified by a software package RNASeq by Expectation Maximization (RSEM). RSEM computes maximum likelihood abundance estimates using the Expectation-Maximization (EM) algorithm for its statistical model, including the modeling of paired-end (PE) and variable-length reads, fragment length distributions, and quality scores, to determine which transcripts are isoforms of the same gene (Li and Dewey, 2011). Fragments per kilobase of exon per million fragments mapped (FPKM) method was used to calculate the expression level as follows: FPKM = [10⁶C/(NL/10³)], where FPKM is the expression of gene A, C is the number of fragments uniquely aligned to gene A, N is the total number of fragments uniquely aligned to all genes, and L is the number of bases on gene A (Ning et al., 2013). The FPKM method can eliminate the influence of different gene lengths and sequencing discrepancy on the calculation of gene expression. The intensity values of each sample were further transformed on log2-scale and used to analyze differential expression. The probe sets with P < 0.01 in at least one of the comparisons were considered as DEGs for further analysis. The WEGO software was used to obtain gene ontology (GO) annotations for the DEGs, while pathway analyses were used to analyze KEGG and BLASTX (E < 0.00001) against the NCBI Nr database (Ye et al., 2006).

The contents of endogenous hormone indole acetic acid (IAA) and gibberellin (GA₃) were determined by enzyme linked immunosorbent assay (ELISA), which was performed as a commercial service in China agricultural university.

Gynophores and pod samples under the same treatments as those in RNA-seq were analyzed using qRT-PCR. cDNAs were reverse transcribed using the PrimeScript™ first-strand cDNA synthesis kit (K1622 Thermo scientific). The PCR was amplified using SYBR Premix Ex Taq™ following the manufacturer's instructions (TaKaRa, Inc., Dalian, China) with the qRT-PCR amplification instrument (ABI 7500, USA). The target gene primers (Table S3) were designed to detect the sample mRNA. Tua5-F and Tua5-R were used as controls to normalize the expression data (Yang et al., 2013). The relative gene expression was calculated using the 2^−ΔΔCT method as described by Livak and Schmittgen (2001). For log₂-transformed FPKM values, the maximum expression level of each selected gene was considered to be 100, and the expression levels of the other genes were transformed accordingly (Zhang et al., 2016).

Ca²⁺ deficiency in soil induces early embryo abortion in peanuts, and this process produces empty pods. This phenomenon was verified by the full pod number and dry weight of pods (Table 1). In our study, we attempted to reveal the molecular mechanism of the effect of free calcium on the development of peanut pods using a RNA-seq approach. A total of 99,030,828 clean reads with an average length of 90 bp, which corresponded to approximately 8.91 gigabase pairs (Gb) of raw data, were obtained. An overview of the sequencing and assembly is outlined in Table S1. 141,819 contigs with an average length of 391 bp were produced. Among these contigs, the length of 81,881 (57.73%) from 100 to 200 bp, 18,847 (13.28%) from 200 to 300 bp, 9256 (6.52%) from 300 to 400 bp, 5633 (3.97%) from 400 to 500 bp, and 26,202 (18.47%) were more than 500 bp in length. The results of the assembly showed 102,819 unigenes with a total length of 102.72 Mb and an average of 999 bp. Of these, the length of 46,924 (45.63%) were from 100 to 500 bp, 17,974 (17.48%) from 500 to 1,000 bp, 12,891 (12.53%) from 1,000 to 1,500 bp, 10,011 (9.73%) from 1,500 to 2,000 bp, and 15,019 (14.60%) were more than 2,000 bp (Figure S1). The length of assembled sequences is a criterion for successful assembly. We calculated the length distribution of contigs and unigenes, and the results showed that the Illumina sequencing solution was reproducible and reliable (Figure S1).

Functional annotation of unigenes gives protein functional annotation, COG functional annotation and Gene Ontology (GO) functional annotation. Unigenes were annotated with the databases of NR, NT, Swiss-Prot, KEGG, COG, and GO. Then the number of annotated unigenes was counted with each database (Figure 1). Analysis of species distribution showed that 42.2% of the unigenes were homologous with the sequence of Glycine max, and 13.8, 13.1 and 11.5% of the unigenes with the sequence of Cicer arietinum, Phaseolus vulgaris, and Medicago truncatula, respectively (Figure 2).

We use Blast2GO program with NR annotation to obtain the GO annotation of the unigenes. GO has three ontologies: molecular function, cellular component and biological process. Every GO-term belongs to a type of ontology (Conesa et al., 2005). Figure S2 showed that the sequences were categorized into 55 functional groups according to the sequence homologies. We mapped the unigenes to the COG database to further predict the possible functions and statistics. In the 25 COG categories, general function prediction (7,524, 32.03%) represented the largest group, followed by transcription (4,202, 17.89%), and replication, recombination, and repair (4,128, 17.57%) (Figure S3).

A total of 14,985 sequences, containing 18,215 potential EST-SSRs, were detected with software MicroSAtellite (MISA). The quantity statistics of SSR classification showed that tri-nucleotide (40.86%) was the most abundant type of repeat motif, followed by di-nucleotide (29.94%), mononucleotide (21.61%), hexa-nucleotide (2.67%), quad-nucleotide (2.53%), and penta-nucleotide (2.38%) repeat units (Table S2). Within the cSSRs detected, AG/CT represented the dominant type, followed by AAG/CTT, ATC/ATG, and AAT/ATT (Figure S4).

Eight DGE libraries were created from four independent biological samples using gynophores and pods with free Ca²⁺-deficient or -sufficient treatments as described in the method section. The clean reads (in millions) obtained from the different libraries were as follows: GD1, 14.91; GD2, 15.59; GS1, 14.54; GS2, 14.34; PD1, 16.02; PD2, 16.67; PS1, 14.19; and PS2, 12.84. The saturation analysis showed that when the number of reads reached a certain amount, the growth curve of detected genes tended to flatten in all libraries, indicating that the sequencing was saturated for gene identification (Figure S5).

GO annotation of the DEGs from two pairwise comparisons (GD/GS and PD/PS) was used to classify genes into different sub-categories belonging to the three main GO categories. In the biological process category, we identified three abundant sub-categories across both comparisons, namely, metabolic process, cellular process, and single organism process. Cell, cell part, and membrane were the main sub-categories identified in the cellular component category, while catalytic activity and binding dominated the molecular function category (Figure 3).

KEGG enrichment analysis allowed the mapping of DEGs to different pathways. The top 20 pathways in each pairwise comparison (mentioned above) are shown in Figure S6. DEGs were mainly enriched in photosynthesis, phagosome, glycolysis/gluconeogenesis, isoflavonoid biosynthesis, glycosylphosphatidylinositol (GPI)-anchor biosynthesis, carbon metabolism, and ABC transporters in GD-VD-GS. Meanwhile, circadian rhythm-plant, plant hormone signal transduction, phagosome, isoflavonoid biosynthesis, glycolysis/gluconeogenesis, fatty acid metabolism, and biosynthesis of unsaturated fatty acids were mainly enriched in PD-VD-PS.

FPKM method was used to calculate the gene expression level. The reference transcripts of a progenitor of cultivated peanut (Arachis ipaensis) and our transcriptome data were used to generate an integrated reference library. For result list of each control-case pair, we draw scatter plots of all expressed genes, different color presents up-regulated, down-regulated or non-regulated genes. The total number of DEGs for each comparison was observed, and a histogram was drawn to show the number of significantly up-regulated or down-regulated genes in each control–case pair (Figure 4).

In order to confirm the transcriptome sequencing results, 10 DEGs in GS, GD, PS, and PD were selected randomly for qRT-PCR to analyze the expression levels. The selected genes participate in calcium signal transduction and phytohormone biosynthesis pathway, material biosynthesis, transcription factor and transporter. The consistency of trends between qRT-PCR and RNA-seq data indicated the credibility of the RNA-seq data (Figure 5).

More than 20 unigenes annotated as transcription factors including NAC, WRKY, bHLH, and ethylene-responsive transcription factors were regulated in GS, while M1/AGAMOUS/DEFICIENS/SRF (MADS) transcription factor family were down-regulated in GS compared with GD. Contrary to the DEGs in pods, no significant difference in calcium-related genes between GD and GS was found, while the expression of some calcium-binding proteins (CBPs), and calcium-binding transcription factors were up-regulated in GS (Table S4).

The expression pattern of lipoxygenases (LOX) genes showed higher expression level in GS and the expression levels of 6 LOX genes were more than 2 times higher in GS than that in GD. Genes related to hormone response were screened out by GO analysis. Several DEGs involved in biosynthesis process and signaling pathway of gibberellin (GA), auxin, ABA, and brassinosteroid were up-regulated in GS at this stage. e. Numerous transporter-related genes were screened from DEGs based on GO analysis. Several energy- and substance-related transporter genes, including sugar transporter, carbohydrate transmembrane transporter, lipid transporter, protein transmembrane transporter, inositol transporter, zinc transporter, sulfate transporter, and nucleobase transporter were all up-regulated in GS. Meanwhile MFS transporter, ABC transporter and myo-inositol transporter had both up-regulated and down-regulated. Only the peptide/histidine transporter was down-regulated in GS (Table S4).

We identified the DEGs between PD and PS considering that calcium plays an important role in the development of peanut pods. Unigenes encoding proteins involved in calcium signaling transduction pathways such as calcium-dependent protein kinases (CDPKs), calcium-binding protein (CBP) and calmodulin (CaM)-binding protein were up-regulated under free Ca²⁺-sufficient conditions. In particular, the expression of Ca²⁺-related protein mitogen-activated protein kinase kinase kinase (MAPKKK) was also higher in PS than that in PD (Table S5).

Two categories of plant hormone-related genes (auxin and GA) showed significantly different expression between PD and PS. The auxin-related genes including 7 auxin response factors (ARFs) and 1 indole-3-acetic acid-amido synthetase were up-regulated in PS. Two selected DEGs, namely, GA 20-oxidase and GA receptor GID1, were up-regulated, while GA 2-oxidase, which is one of the key enzymes that catalyzes the inactivation of biological GA, was down-regulated in PS. Meanwhile, the ent-kaurenoic acid oxidase (KAO), an important enzyme in GA biosynthesis (Yamaguchi and Kamiya, 2000; Paparelli et al., 2013), was up-regulated under free Ca²⁺-sufficient conditions. These results were in agreement with the high contents of IAA and GA₃ in PS (Figure 6). Thus, our results indicated that the up-regulated calcium-related genes increased GA level, which may play an important role in the early embryo development.

Genes encoding storage protein oleosin and some synthetases, such as tetrahydrocannabinolic acid synthase, long chain acyl-CoA synthetase, and callose synthase, exhibited higher level (5–8-folds) under free-Ca²⁺-sufficient conditions. In addition, eight lipoxygenase were also up-regulated in PS. Previous studies indicated that lipoxygenase is the first key enzyme in the synthesis pathway of jasmonic acid, which plays an important role in the germination, growth, and stress resistance of plant seeds (Rahimi et al., 2016). GIGANTEA, a gene involved in circadian clock and phytochrome signaling (Cha et al., 2017), was down-regulated in PS (Table S5). This phenomenon may be due to the fact that this light response gene cannot play a regulatory role in the subterranean areas, because darkness is also a prerequisite for the normal development of peanut pods.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.