Peanut cultivar, ‘Huayu 22,’ was provided by Biotechnology Research Center, Shandong Academy of Agricultural Science (SAAS, China). Free Ca²⁺ content in soil is 14 g (free Ca²⁺) kg^-1 (dry soil). The pods were collected at the 1st (Stage 1, abbr. S1), 5th (Stage 2, abbr. S2), 10th (Stage 3, abbr. S3), 20th (Stage 4, abbr. S4) day after the peg elongation into the soil. Four stage materials were used for Ca²⁺ excretion determination. The materials of S2, S3, and S4 were collected and immediately frozen in liquid nitrogen, and then stored in a freezer at -80°C for comparative transcriptome analysis.

Field experiments was conducted to observe effects of free Ca²⁺ in soil on pod development in Sanzhuang, Rizhao, Shandong Province in spring 2014 and 2015. Free Ca²⁺ content in soil is 4 g (free Ca²⁺) kg^-1 (dry soil). The 120 kg CaO ha^-1 was applied in soil as Ca²⁺ treatment, and non-CaO treatment was as control. Each treatment had three replicates. CaO was applied before sowing.

Ca²⁺ excretion in pods of peanut were measured with NMT system (NMT100 Series, YoungerUSA LLC, Amherst, MA, United States; Xuyue (Beijing) Science and Technology Co., Ltd., Beijing, China) and iFluxes/imFluxes 1.0 (YoungerUSA, LLC, Amherst, MA, United States) Software. Steady-state ion fluxes were measured for 5–20 min. After that, the test treatment was applied and the ion flux in the meristematic zone (120 μm from the root/pod tip) was measured for further 10 min and recorded the data (Yue et al., 2012). And the roots and pods were treated with 40 μm Ca²⁺ solution. Each treatment was repeated at least three times and all tests were repeated at least six times.

Total RNA was isolated and integrity confirmed using a 2100 Bioanalyzer (Agilent Technologies). Beads with oligo(dT) were used to isolate poly(A) mRNA from total RNA (Qiagen GmbH, Hilden, Germany).

Following purification, the mRNA was fragmented using divalent cations under elevated temperature. Taking these short fragments as templates, the first-strand cDNA was synthesized using random hexamer primers and Superscript^TM III (Invitrogen^TM, Carlsbad, CA, United States). The second strand cDNA was synthesized using buffer, dNTPs, RNaseH, and DNA polymerase I. Short fragments were purified with a QiaQuick PCR extraction kit (Qiagen) and resolved with EB buffer for end reparation and poly(A) addition. The short fragments were then connected using sequencing adapters. After agarose gel electrophoresis, suitable fragments were used as templates for PCR amplification. Finally, the library was sequenced from both directions on HiSeq 2000 System (illumina, San Diego, CA, United States) with 100 bp of data collected per run by applying TruSeq PE Cluster and TruSeq SBS Kits (illumina, San Diego, CA, United States). Data analysis and base calling were achieved by applying the illumina instrument software.

The reference transcripts of a progenitor of cultivated peanut (Arachis ipaensis) were used to generate an integrated reference library. The processing and assembly of transcriptome data were performed with the application of modified Velvet to construct unique consensus sequences (Zerbino and Birney, 2008). The gene expression profile was developed by mapping trimmed transcriptome reads onto the unique consensus sequences by using SOAP2 (Li et al., 2009).

Comparative transcriptomics analysis was conducted among the three different developmental stages. The intensity values of each sample were further transformed on log₂-scale and used for differential expression analysis. The probe sets with a P-value < 0.01 and > two-fold changes in at least one of the comparisons were considered as differentially expressed genes (DEGs) for further analysis. Unigenes were used for BLAST searches and annotation against an NCBI Nr protein database (NCBI non-redundant sequence database) by using an E-value cut-off of 10^-5 (E-value < 0.00001). Unigenes sequences were further aligned by BLASTX to protein databases such as Swiss-Prot, KEGG, and COG, retrieving proteins with the highest sequence similarity with the given unigenes along with their protein functional annotations. If results of different databases conflicted, a priority order of Nr, Swiss-Prot, KEGG, and COG was followed. For unigenes that did not align to any of the above databases, EST Scan software (Iseli et al., 1999) was used to predict their coding regions and determine sequence direction. Unigenes aligned to databases with higher priority did not enter the next circle. The alignments were considered complete when all circles were finished. The coding region sequences were then determined for proteins with the highest ranks using BLAST. Unigenes that could not be aligned to any database were scanned by EST Scan (Iseli et al., 1999) to determine the nucleotide (50–30) and amino acid sequences of the coding regions. The Blast2 GO was used to obtain Gene Ontology (GO) annotations for the unigenes, as well as for KEGG and COG analysis (Conesa et al., 2005). The WEGO software was then used to perform GO functional classification of all unigenes to view the distribution of gene functions of the species at the macro level (Ye et al., 2006). It was mapped all of the annotated unigenes to GO terms in the database and calculated the number of unigenes associated with every term.

Total RNA was isolated and purified using TRIzol Reagent (Tiangen, China) according to the manufacture’s instruction. The first-strand cDNA was synthesized by using Reverse Transcriptase M-MLV (Takara, Japan) according to the manufacturer’s protocol. The PCR was amplified following the instruction of SYBR Premix Ex Taq^TM (Takara, Japan) with the fluorescent quantitative PCR amplification instrument (ABI 7500, United States). The target gene primers were used to detect the sample mRNA. Control reactions were carried out by using primers Tua5-F and Tua5-R to detect the transcript encoding ubiquitously expressed tua5 (Chi et al., 2012). All assays for a target gene were performed in triplicate synchronously under identical conditions.

Peanut is one of the important oil crops for China, and its yield is often affected by all kinds of environmental factors, including free Ca²⁺ content in soil. When free Ca²⁺ is deficient in soil, it will cause unfilling of pods. Figure 1 showed that pods were unfilled under 4 g (free Ca²⁺) kg^-1 (dry soil), and CaO application can benefit pod filling.

To assess the role of Ca²⁺ on pod development, hair root and four stages of pods were analyzed to determine Ca²⁺ excretion where: the peg just entered the soil (S1), tip of the peg initiated expansion (S2), beginning of pod reticulation (S3), and when the pod was clearly reticulated (S4), respectively, used to analyze (Figure 2A). The non-invasive microtest technique was used to detect net fluxes of Ca²⁺ at different seed tissues of different stages. Besides hair roots, peanut shells showed Ca²⁺ excretion at all four stages; however, both the embryo of S4 and the red skin of S3 and S4 showed Ca²⁺ absorbance. These results showed that embryo alone and embryo with red skin of peanut need Ca²⁺ during development (Figure 2B).

To further detect some mechanisms relate to the roles of Ca²⁺ in pod development. About 70 million high quality reads were obtained from the library. The average read size, Q30 percentage (sequencing error rate, 0.1%), and GC percentage for each library was 90 bp, >95%, and >40%, respectively. Clean reads from each library were used for assembly separately. And 263,110 unigenes were generated from library. The de novo assembly of all sequencing data using the modified Velvet was followed by the application of the SOAP2 program for unigene identification. The average sequence length of the unigenes was 1,235 bp and its range was 201–4,000 bp (Figure 3A).

The GO assignments were used to categorize the functions of the predicted peanut Unigenes. Based on the sequence homologies, the sequences were categorized into 45 functional groups (Figure 3B), contained three main categories of the GO classification (biological processes, cellular components, and molecular functions). The genes involved in cellular process and metabolic process were dominant in the “Biological Process” category. “Cell,” “Cell part,” and the “Organelle” are the top three abundant categories in “Cellular component,” While “Binding,” and “catalytic activity” are dominant in the “Molecular function” category (Figure 3B).

To investigate expression levels of stage-specific genes during pods development, we conducted comparative analysis of the transcriptome profiles among three pod development stages at 5th (S2), 10th (S3), 20th (S4) day after the peg elongation into the soil. As shown in Figure 4A, totally 1,169 DEGs were obtained between S2 and S3, including 50 up-regulated and 1,119 down-regulated genes. Similarly, 1,172 DEGs (including 53 up-regulated and 1,119 down-regulated genes), 7 DEGs (including two up-regulated and five down-regulated genes) were identified between S3 and S4, S2 and S4, respectively. In addition, some other stage-specific expression genes were also analyzed. Approximately 39, 2157, and 1 stage-specific expression genes were identified between S2 and S3, S3 and S4, S2 and S4, respectively (Figure 4B).

Since deficient free Ca²⁺ in soil can cause unfilling of pods, some DEGs involved in biological process of Ca²⁺ binding and Ca²⁺ signaling were identified based on GO analysis. These potential candidate genes contained 53 Ca²⁺ or Ca²⁺ signal pathway related genes (Table 1). These genes were significant differentially expressed during seed and pod development, suggesting that Ca²⁺ and calcium related protein might be involved in seed development.

Genes related to hormone response were screened out by GO analysis. There were 40 auxin related genes, 15 gibberellin related genes, 20 ethylene related genes, 2 abscisic acid related genes, and 7 cytokinin related genes (Table 2) showed significant differentially expression during seed and pod development.

In order to confirm the transcriptome sequencing results, Ca²⁺ and hormone related genes were selected from DEGs based on the GO analysis and used to real-time fluorescent quantitative PCR analysis. As shown in Figure 5, the expression pattern of 8 selected DEGs, which are related to Ca²⁺ signal pathway, were consistent with their respective microarrays data. Among these genes, two genes, which includes one CDPK gene (comp33153) and one uncharacterized calcium-binding protein At1g02270 (comp46176), were up-regulated at S4. The Serine/threonine-protein kinase SRK2E (comp45239) was up-regulated at S2, down-regulated at S3, and then up-regulated at S4. Three probable calmodulin like protein CMLs (comp33933, comp43665, and comp50658), one calcineurin B-like protein (comp47042) and one calcium-binding protein PBP1 (comp43681) were up-regulated at S2 and then down-regulated during the later stages.

Among the 12 selected auxin related genes, six genes, i.e., auxin-induced protein 5NG4 (comp43478), auxin-responsive protein IAA31 (comp35718), auxin-induced in root cultures protein 12 (comp50484), auxin-induced protein 15A (comp35357), auxin-induced protein 10A5 (comp36019), and auxin transporter-like protein 3 (comp36785), were down-regulated in S2 and then up-regulated during the later stages. Three auxin response factor 8 (comp52994, comp45289, and comp46658), one auxin-responsive protein IAA26 (comp36524) were down-regulated at S4. One auxin-induced protein 22E (comp45084) was down-regulated at S2, up-regulated at S3, and then down-regulated at S4 (Figure 6).

Among 11 selected gibberellin related genes, four genes, i.e., two gibberellin-regulated protein (comp41469, comp27867), two gibberellin 2-beta-dioxygenase (comp45036, comp43715), were down-regulated at S2 and then up-regulated at S4. One transcription factor MYB44 (comp40799), one gibberellin 20 oxidase (comp44953) and one gibberellin 3-beta-dioxygenase (comp44039) were up-regulated at S4. One gibberellin 20 oxidases (comp53643), one gibberellin 2-beta-dioxygenase (comp47582), one gibberellin 3-beta-dioxygenase (comp45509) and one gibberellin-regulated protein (comp13386) were up-regulated at S2, down-regulated at S3 and then up-regulated at S4 (Figure 7).

Seventeen selected ethylene related genes were analyzed. Two ethylene-responsive transcription factor RAP2 (comp53691 and comp42029), one ethylene-responsive transcription factor WRI1 (comp48452), three ethylene-responsive transcription factor ERF (comp34495, comp37722, and comp11865) and one AP2-like ethylene-responsive transcription factor AIL (comp47962) were down-regulated at S2 then up-regulated at S4. One ethylene-responsive transcription factor ERF (comp35300), one ethylene-responsive transcription factor (comp44242), two ethylene-responsive transcription factor CRF (comp47316 and comp43675) and one AP2-like ethylene-responsive transcription factor ANT (comp43442) were only down-regulated at S4. The AP2-like ethylene-responsive transcription factor TOE (comp47833) and ethylene-responsive transcription factor ERF (comp45116) were up-regulated at S4. AP2-like ethylene-responsive transcription factor At2g41710 (comp42987) and AP2-like ethylene-responsive transcription factor AIL1 (comp48832) were down-regulated at S2, up-regulated at S3, and then down-regulated at S4. The ethylene-responsive transcription factor 13 (comp37530) was up-regulated at S2 and S3, then down-regulated at S4 (Figure 8).

Among the four cytokinin related genes, cytokinin dehydrogenase (comp42347) was only up regulated at S4, cytokinin hydroxylase (comp44683) was down-regulated at S2 then up-regulated in later stage, cytokinin riboside 5′-monophosphate phosphoribohydrolase LOG3 (comp53620) was down-regulated from S3 and the Cytokinin riboside 5′-monophosphate phosphoribohydrolase LOG8 (comp38806) up-regulated from S3 (Figure 9A).

For the two abscisic acid related genes, the abscisic acid 8′-hydroxylase 3 (comp29840) was up-regulated at S3 and then down-regulated, however, the abscisic acid insensitive 5-like protein 1 (comp38442) was down-regulated at S2 (Figure 9B).