Transcriptome analysis reveals the mechanism of improving erect-plant-type peanut yield by single-seeding precision sowing

Plant materials and growth conditions

Peanut (Arachis hypogaea L.) cultivar “Huayu 22” (provided by the Shandong Peanut Research Institute, China), a large-grain peanut cultivar, was used in this study. Peanut seeds were cultivated under field pot conditions at the Shandong Academy of Agricultural Sciences Station (117°5′ E, 36°43′ N), Ji’nan, China. The soil is sandy loam, containing 1.1% organic matter (W/W), 82.7 mg kg⁻¹ alkali-hydrolyzed nitrogen, 36.2 mg kg⁻¹ available phosphorus, 94.5 mg kg⁻¹ available potassium and 14.9 g kg⁻¹ exchangeable calcium. The SS sowing treatment consisted of 237,000 holes per hectare with 1 grain per hole, while the DS sowing consisted of 138,000 holes per hectare with two grains per hole as the control; the distances between two adjacent holes were 10.5 cm and 18 cm in the SS and DS treatments, respectively. Three biological replicates were used in this study. The mixed tissues were used in transcriptome sequencing.

RNA extraction and cDNA library construction

When the peanut plants reached the seedling stage, the top third leaf and the whole root systems were collected for RNA extraction and library construction. Total RNA was isolated from the leaves and roots respectively, using the total RNAiso Reagent (TaKaRa, Dalian, China) according to the manufacturer’s instructions. The quality and purity of RNA samples were detected using the Agilent 2100 Bioanalyzer platform and Agilent RNA 6000 Nano Kit (Agilent, Santa Clara, CA, USA). Qualified RNA samples were used to isolate mRNA with the oligo (dT) method and the mRNAs were fragmented. Then, first-strand cDNA and second strand cDNA were synthesized. Next, the purified cDNA fragments were linked with adapters, and sequencing was performed by a commercial service provider (BGI Tech, Shenzhen, China). The sequenced reads that contained adaptor sequences, low-quality bases, and high contents of unknown bases (i.e., N calls) were removed before downstream analyses. After read filtering, clean reads were mapped to the Arachis_duranensis reference genome (BioSample: SAMN02982871, BioProject: PRJNA258023) using Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT). The sequence data can be located from https://www.ncbi.nlm.nih.gov/assembly/GCF_000817695.2/. Then, Genome Analysis Toolkit (GATK) was used to call single SNPs and INDELs (McKenna et al., 2010) and RNA-Seq by Expectation-Maximization (RSEM) was used to calculate gene expression levels for each sample, which were then determined using the fragments per kb per million (FPKM) mapped fragments method developed by Li & Dewey (2011). DEGseq algorithms were used to assess the DEGs. Based on the gene expression level, DEG was identified between samples or groups. Gene Ontology (GO) was performed classification and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway classification based on these DEGs. The GO framework includes three ontologies: molecular biological function, cellular component and biological process. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. Co-expression analysis was performed using the WGCNA R package according to the methods detailed by previously published studies (Gao et al., 2018; Song et al., 2018).

Yield composition per plant and measurement

At maturity, 10 representative plants were selected to investigate the number of pods, full pods and double kernels, as well as pod weight per plant. The pods and plants were dried to constant weights, and the economic coefficient was calculated as follows: economic coefficient = pod dry weight/(plant dry weight + pod dry weight).

Antioxidant enzyme activities and root activity

Leaf tissues (0.5 g) were ground up with phosphate buffer (pH 7.8) containing 0.1 mM EDTA and 1% (g ml⁻¹) PVP and centrifuged at 4 °C for 10 min. Then, the supernatants were used as enzyme extracts. The nitroblue tetrazole (NBT) method was used to detect SOD activity. The absorbance value was determined at a 560 mm wavelength after the reaction; 50% inhibition of NBT reduction was regarded as the enzyme activity unit (U), expressed as U g⁻¹ FW. The guaiacol method was used to determine the POD activity; the increase in OD₄₇₀ per min was used as the enzyme activity unit (U), expressed as Δ470 g⁻¹. CAT activity was also determined by the guaiacol method. The decreased in OD₂₄₀ per min was used as the unit of enzyme activity (U) and expressed as mg·g⁻¹·min⁻¹ (Giannopolitis & Ries, 1977; Aebi, 1984; Jimenez et al., 1997). Root activity was measured by the TTC method; 0.5 g of root tip samples were revolved in a mixture of 0.4% TTC and phosphoric acid buffer solution and kept in the dark for 1–3 h at 37 °C. After that, 2 ml of 1 mol L⁻¹ sulfuric acid was added to stop the reaction. Simultaneously, a blank experiment was conducted. The root samples and sulfuric acid were added first, and 10 min later, the other reagents were added. The operation was the same as that described above. The roots were removed, dried, and ground together with 3–4 ml of quartz sand and ethyl acetate in a mortar. The supernatants were transferred to a test tube, and the residue was washed with a small amount of ethyl acetate two or three times, with all ethyl acetate collected into the test tube. Finally, ethyl acetate was added to bring the total volume up to 10 ml. Blank experiments were used as the control, and the OD value of the reaction solution was measured at λ = 485 nm.

qRT-PCR verification analysis

Leaves and roots of plants grown under the SS and DS treatments, as in the RNA-seq experiment, were selected randomly and analyzed using qRT-PCR. The qRT-PCR amplification instrument (ABI 7500 fast; Applied Biosystems, Foster City, CA, USA) was used to amplify the related genes with SYBR Premix Ex Taq^™ (TaKaRa, Nugegoda, Sri Lanka) following the manufacturer’s instructions (Yang et al., 2013). The relative gene expression was calculated using the 2^−ΔΔCt method described by Livak & Schmittgen (2001). For log₂-transformed FPKM values, each selected gene’s maximum expression level was considered to be 100, and the expression levels of the other genes were transformed accordingly (Zhang et al., 2016).

Determination of Res content

The Res extraction processes were carried out according to the method detailed by Tang et al. (2010) with minor modification; 0.25 g of leaves or roots were collected and ground into a fine powder, which was dissolved with 20 ml of 95% ethanol in a flask and oscillated at room temperature overnight. After centrifugation at 8,000 rpm for 10 min, the supernatants were subjected to HPLC analysis. The extracts were filtered through a 0.45-μm membrane and then tested with the Rigol L3000 HPLC system (Rigol Technologies, Beaverton, OR, USA) using a Kromasil C18 reversed phase column (4 μm, 250 mm × 4.6 mm; Kromasil, Bohus, Sweden) at room temperature. Mobile phase A consisted of 0.1% phosphoric acid–H₂O, while phase B consisted of acetonitrile (A:B = 70:30, flow rate 1 ml min⁻¹); the determination wavelength was 280 nm.

Transcriptome sequencing and gene analysis

The peanut plants grown under different planting patterns differed little in their above-ground portions, while great differences were observed between their root systems. The basis of high and stable crop yields is high biological yields, which depends largely on the development of root systems. The capillary lateral root is the main type of root exhibited under SS. Thus, the capillary lateral root length and its proportion within the root system were promoted accordingly, thereby increasing the total root length and root absorption area of peanut plants. This phenomenon was verified by the mean absorbing area and dry weight of root tissue observed under SS (Fig. 1).

A total of 12 cDNA libraries were constructed from leaf and root tissues of peanut plants grown under SS and DS sowing treatments. All 12 samples were sequenced on the Illumina HiSeq Platform in total, generating about 6.55 Gb per sample. The reads that contained adaptor sequences of low-quality or had high contents of unknown base (i.e., N) calls were removed, resulting in a total of 523,800,338 clean reads that were acquired, with an average of 43 million reads per libraries (Table 1). Approximately 77% and 80% of reads from leaf and root tissues were mapped to the reference genome and 36,778 genes were identified, of which 34,529 are known genes and 2,322 are novel genes. The RNA-seq sequencing data for the present work has been uploaded in NCBI Sequence Read Archive under BioProjects, PRJNA497502 (SRA: SRP166140).

Single nucleotide polymorphisms (SNP), which include transitions or transversions of single bases, refer to the difference of a single nucleotide (A, T, C, or G) between homologous DNA sequences, a key element of the diversity of genomes among species or individuals. There were more transitions than transversions among all the samples. Among SNPs, transition, transversion, A–G and C–T were the most abundant terms (Fig. S1). Insertion–Deletion (INDEL) refers to the insertion or deletion of the small fragments (one or more, less than 50 bp) that occur in a sample relative to the reference genome. Most INDELs occurred in exonic and intronic sequences, with the proportion differing among samples (Fig. S2). Based on the results of SNP, INDEL and gene expression, our study was presented in the form of a ring diagram with Circos software (Fig. S3).

Differentially expressed gene detection

A total of 2,567 and 2,706 differentially expressed genes (DEGs) were identified from leaves to roots, respectively. Under SS sowing treatment, 544 and 1,771 genes were expressed at a higher level, while 2,023 and 935 genes were expressed at lower levels compared with the DS sowing treatment in leaves and roots, respectively (Fig. 2). The expression patterns of the majority of DEGs differed between roots and leaves. For example, only 83 (15.2%) of the 1,771 genes up-regulated in roots were also up-regulated in leaves (Fig. 3). Notably, some genes in roots and leaves even showed opposite expression patterns; for example, forty DEGs were up-regulated in leaves but down-regulated in roots. However, 245 DEGs in roots and leaves had the same expression trends, including 83 up-regulated DEGs and 162 down-regulated DEGs (Fig. 3).

qRT-PCR verification of RNA-seq results

To verify the RNA-seq data, quantitative real-time PCR (qRT-PCR) was used to test the expression of 24 genes with different functional assignments (Figs. 4A and 4B). Among them, six and five genes were up-regulated under the SS sowing treatment in leaves and roots, respectively. UDP-glycosyltransferase, chalcone synthase and GPI mannosyltransferase were each encoded by one of these genes, while two genes encoding hypothetical proteins were also used for qRT-PCR detection. The other eleven selected genes were downregulated both in roots and leaves. These genes included one phospholipase, one seed linoleate 9S-lipoxygenase, one zinc finger SWIM domain-containing protein, and one pectinesterase. The RNA-seq data were consistent with the qRT-PCR results for these 24 genes. The correlation between the relative expression (log₂ SS/DS) estimated by RNA-seq and the qRT-PCR results were rather high in leaves (R² = 0.9772) and a little lower in roots (R² = 0.8884; Figs. 4C and 4D).

Functional analysis of DEGs

Blast2GO was used to perform GO classification and functional enrichment (Conesa et al., 2005) on the identified DEGs. The sequences were categorized into 45 functional groups according to sequence homology. The main categories of biological process, molecular function, and cellular component were visualized employing WEGO (Fig. 5). In the biological process category, “metabolic process”, “single-organism process” and “cellular process” terms were enriched, suggesting the metabolic activity is higher under SS sowing. In the cellular component category, “cell”, “cell part” and “membrane” were the most abundant terms. The most enriched category was molecular function, with enrichment for “binding”, “catalytic activity” and “transporter activity” in particular, suggesting a high level of metabolic activity changes under SS sowing (Fig. 5).

Kyoto Encyclopedia of Genes and Genomes pathway classification was performed to acquire biological information for understanding the regulatory networks and molecular mechanisms associated with the SS sowing treatment. DEGs were mainly enriched in the MAPK signaling pathway, glycerolipid metabolism, phenylalanine metabolism, sphingolipid metabolism, isoflavonoid biosynthesis, flavonoid biosynthesis and tryptophan metabolism in leaves. Meanwhile, phenylpropanoid biosynthesis, biosynthesis of secondary metabolites, MAPK signaling pathway, flavonoid biosynthesis, zeatin biosynthesis and flavone and flavonol biosynthesis were mainly enriched among DEGs in roots (Fig. S4). Interestingly, all of these pathways participated in biosynthesis associated with particular metabolic processes, suggesting that these processes were activated. Weighted gene correlation network analysis (WGCNA) is a systematic biological method for describing gene association pattern among different samples. Information about nearly 10,000 genes, corresponding to the genes with the greatest changes in expression, was used to identify gene sets of interest and analyze associations with phenotypes. In our study, the genes from all the samples were divided into modules for analysis, with twelve colors representing each of the different modules in Fig. 6.

Changes of Res synthesis related pathway between cultivation techniques

In our study, a total of 20 Res synthesis-related genes were dramatically induced, while no Res synthesis-related genes were down-regulated in the roots under the SS sowing treatment. In contrast, a total of 10 DEGs related to Res synthesis were down-regulated in the leaves (Table 2). We also measured the Res content, and the results were in agreement with the gene expression levels in leaf and root tissues (Fig. 7). Accordingly, compared with the roots of plants grown under DS sowing, the higher Res content under SS precision sowing appeared to support the initiation of the defense response against toxins produced by Aspergillus flavus; pathogenesis-related proteins, TFs (e.g., WRKY, bZIP, ERF) and 4CL were up-regulated (Table 3).

Key genes related to stress tolerance under SS

Transcription factors, which bind to the cis-elements upstream of promoters, have been reported to orchestrate abiotic stress responses (Joshi et al., 2016). Under the SS precision sowing treatment, the higher expression of these TF genes (WRKY, MYB, bZIP, ERF) may improve stress tolerance compared with that observed in the DS treatment. In plants, redox processes also play an important role in stress tolerance. Several oxidoreductase genes encoding ascorbate oxidase, glutaredoxin, and cytochrome P450 monooxygenases (CYP) exhibited increased expression in roots under SS relative to DS. Specifically, compared with DS, about thirty genes encoding CYP were up-regulated under SS (Table 3). In addition, four genes encoding l-ascorbate oxidase homologs, which take part in ascorbate recycling, were enriched in the SS sowing treatment. In plants, ascorbate contributes to improving tolerance against various stresses by regulating the levels of cellular H₂O₂ (Ishikawa & Shigeoka, 2008). Higher enzyme activities of POD, SOD, and CAT were also observed in roots under SS compared to that under DS (Table S1). Four genes encoding glutaredoxin, which might promote reducing disulfide bridges, were slightly induced under SS compared with DS (Table 3).

Peanut pod yields and associated genes under the SS sowing pattern

The number of pods and pod weight per plant is the direct factors affecting peanut yield. In this study, yield estimation indicated that the number of pods, full pods, and double kernels per plant under SS was higher than under DS. The theoretical pod yield of SS reached 1.225 kg, which was 12.9% higher than that of DS (1.085) (Table 4). Hormones, such as indole-3-acetic acid (IAA), abscisic acid (ABA), gibberellin (GA) and brassinosteroid (BR), play important roles in the regulation of growth and development in plants (Tian et al., 2017). Auxin response factor (ARF), a TF involved in auxin and regulating plant root growth and seed development (Salmon, Ramos & Callis, 2008), showed increased expression levels under SS. Nine genes that participate in auxin and BR biosynthesis were abundantly expressed under SS compared with under DS (Table 3). Two lipoxygenase family genes also had increased expression levels under SS.