Transcriptome profiling provides insights into molecular mechanism in Peanut semi-dwarf mutant

Background Plant height, mainly decided by main stem height, is the major agronomic trait and closely correlated to crop yield. A number of studies had been conducted on model plants and crops to understand the molecular and genetic basis of plant height. However, little is known on the molecular mechanisms of peanut main stem height. Results In this study, a semi-dwarf peanut mutant was identified from 60Co γ-ray induced mutant population and designated as semi-dwarf mutant 2 (sdm2). The height of sdm2 was only 59.3% of its wild line Fenghua 1 (FH1) at the mature stage. The sdm2 has less internode number and short internode length to compare with FH1. Gene expression profiles of stem and leaf from both sdm2 and FH1 were analyzed using high throughput RNA sequencing. The differentially expressed genes (DEGs) were involved in hormone biosynthesis and signaling pathways, cell wall synthetic and metabolic pathways. BR, GA and IAA biosynthesis and signal transduction pathways were significantly enriched. The expression of several genes in BR biosynthesis and signaling were found to be significantly down-regulated in sdm2 as compared to FH1. Many transcription factors encoding genes were identified as DEGs. Conclusions A large number of genes were found differentially expressed between sdm2 and FH1. These results provide useful information for uncovering the molecular mechanism regulating peanut stem height. It could facilitate identification of causal genes for breeding peanut varieties with semi-dwarf phenotype.

developing semi-dwarf peanut variety are of great significance.
Stem elongation is determined by cell growth including cell division and cell expansion. The dynamic variation of plant cell wall is critical in these processes [3]. The cellulose synthase super-family containing cellulose synthases (CesA) and cellulose synthase-like genes (CSL) is very important in cell-wall biogenesis and remodeling. CesA or CSL mutation in plants resulted in a dwarf phenotype with short internode or hypocotyl and small leaves [4][5][6][7][8]. Cell wall extensibility controls the rate of plant cell growth, especially cell expansion [9]. Expansins are plant cell-wall loosening proteins which promote cell enlargement by influencing cell wall extensibility and are essential in many critical developmental processes [10]. Suppression of expansin genes resulted in small, curled leaves, reduced plant height and early flowering [11][12][13][14]. Besides expansins, xyloglucanendotransglucosylase/ endohydrolases (XTHs) were also suggested to regulated cell wall loosening [3,15]. XTHs catalyze the endo cleavage of xyloglucan polymers and subsequent restructuring, which allow cellulose microfibrils moving along with the cell expansion or elongation process [16][17][18]. EGases (endo-1,4-β-glucanases) have been proposed to modify hemicellulose network, and were involved in processes that require cell wall weakening, including cell elongation, organ abscission, and fruit softening [19].
Auxin polar transport is critical in regulating the whole progress of plant growth and development [49]. AUX1/LAX family, PIN family and ABCB subfamily are three classes of polar auxin transporters. MDR (Multidrug Resistance) genes encode P-glycoproteins and belong to subfamily B of ATP-binding cassette (ABCB). The atmdr1 mutants were defective in basipetal auxin transport and shorter than wild-type plants [50]. SAUR (Small Auxin up RNA) genes, the largest family of auxin-responsive genes, functioned as positive effectors of cell expansion and were up-regulated in auxinmediated cell elongation [51][52][53].
Many transcription factors are important regulators of plant height. SHORT INTERNODES (SHI) gene is a negative regulator of GA response and its overexpression mutant shi showed a dwarf phenotype similar to mutants defective in GA biosynthesis or response [54]. WRKY46/54/70 could be activated by BRs and were cofactors of BES1 to modulate BR target genes [55]. PACLOBUTRAZOL RESISTANT (PRE) is a positive regulator of cell elongation which can be activated by BRs, GAs and repressed by light [56,57]. PRE together with IBH1, HBI1 and ACE constitute a triantagonistic cascade module which integrate hormonal and environmental signals to regulate cell elongation [58][59][60].
In peanut, genetic basis underlying plant height remains unclear. Some studies to identify QTLs for plant height have been conducted in peanut. Huang [63]. However, no gene closely related to plant height has been cloned in peanut until now.
In the previous research, we built a mutant library through 60 Co radiation. Among which, two mutants displayed semi-dwarf phenotype. In this study, genomewide gene expression profiling of stem and leaf from sdm2 was conducted to understand the underlying regulation network of main stem height in peanut. The expression of several genes involved in biosynthesis and signal transduction of BRs and GAs were drastically altered in sdm2. Accordingly, the downstream genes related to cell wall biogenesis and remodeling were also changed at the transcriptional level.

Phenotype characteristics of peanut dwarf mutant
The main stem height of sdm2 is significantly shorter than that of FH1 (P < 0.01) (Fig. 1a). We measured the main stem height of mutant and FH1 grown in the experimental farm every 2 weeks during the vegetative growing stage. Results showed that during the first 2 weeks, main stem height between sdm2 and FH1 had no obviously difference. At 28 DAP (days after planting), main stem height of sdm2 was significantly shorter than that of FH1 (P < 0.01). In the following developmental stages, the height difference between sdm2 and FH1 became more and more obvious. At 84 DAP, the sdm2 was 15.9 cm in height, only about 60% of FH1 (26.8 cm) (Fig.  1b). The leaf of sdm2 was smaller, thicker and dark green in color, and showed delayed senescence compare to that of FH1. At harvest time, the number of internodes of sdm2 was 3~4 less than FH1. The number of branches (P < 0.05) and pods (P > 0.05) per plant was obviously more, while the 100-pod weight and 100-kernel weight were significantly decreased in sdm2 (P < 0.01).
To further accurately measure the change of mutant main stem height, we conducted hydroponic culture of sdm2 and FH1 in Hogland solution. The length of each internode and hypocotyl of 14 DAP seedlings were measured. We found that it was different from the results in experimental farm. For 14 DAP seedlings, the sdm2 height was already significantly shorter than that of FH1 (P < 0.05) (Fig. 1c, d). There were total five internodes in FH1 and four in sdm2 seedlings. The hypocotyl length of sdm2 was slightly shorter than FH1 and there was no significant difference (P > 0.05) for the 14 DAP seedlings. The length of most internodes of sdm2 was significantly shorter than FH1 (Fig. 1d). These results indicated that the height difference between sdm2 and FH1 was caused by both reduced total number of internode and length of each internode.

Sequencing data analysis
To clarify global gene expression changes in sdm2, 12 cDNA libraries were constructed with stem and leaf from mutant and FH1 seedlings. The libraries were sequenced using BGISEQ-500 platform. A total of 0.26 billion raw reads were generated, and the average output of each sample was 21.7 million. After removing adaptor sequences, low-quality and N-containing reads, an average of 21.58 million clean reads were obtained from each library with the clean read ratio of 99.47% (Table 1). Approximately 19.27 million and 16.00 million clean reads in each library matched the reference genome and gene set perfectly with the average mapping ratio of 89.31 and 74.18%, respectively. To investigate gene expression correlation among samples, the Pearson correlation coefficient of all gene expression levels between each two samples were calculated (Additional file 1: Fig. S1). It showed a higher correlation of similar tissue between sdm2 and FH1. The principal component analysis (PCA) could clearly divide all samples into two clusters according to two kinds of tissues (Additional file 2: Fig. S2).

DEGs between sdm2 and FH1
From all the samples, a total of 65,435 genes were identified. The gene expression of stem and leaf between sdm2 and FH1 was analyzed, and 3733 and 3715 differentially expressed genes (DEGs) were identified, respectively (Additional files 3 and 4: Table S1 and S2). In sdm2 stem, 1786 genes were up-regulated and 1947 genes were down-regulated to compare with FH1 ( Fig. 2a and Additional file 3: Table S1). In sdm2 leaf, the up-and down-regulated genes were 1644 and 2071, respectively ( Fig. 2a and Additional file 4: Table S2). The majority of DEGs had different expression pattern in stem and leaf. Among 3733 DEGs in stem, only 873 were also differentially expressed in leaf. There were 473 down-regulated and 163 up-regulated DEGs with the same change trend in stem and leaf. It is worth noting that the number of common down-regulated DEGs was obviously higher than that of up-regulated ones. Interestingly, a large number of DEGs showed opposite expression pattern in stem and leaf. For example, there were 149 genes downregulated in mutant stem while up-regulated in leaf and 88 genes up-regulated in mutant stem but downregulated in leaf (Fig. 2b).

Functional analysis of DEGs
GO classification showed that 1724 DEGs from stem and 1716 DEGs from leaf were classified into three categories: biological process, cellular component, and molecular function ( Fig. 3 and Additional file 5: Table S3). In molecular function category, catalytic activity and binding were the most abundant terms. For cellular component  category, cell, membrane, membrane part and organelle were the main terms. The top two terms of biological process were metabolic process and cellular process. The top 20 terms in stem and leaf were shown in Additional file 6: Fig. S3 and Additional file 7: Table S4. KEGG pathway analysis was also carried out. According to the enrichment results, the top 20 pathways in stem and leaf were shown in Fig. 4 and Additional file 8: Table S5. In stem, many pathways involved in hormone biosynthesis and signal transduction including terpenoid backbone biosynthesis, monoterpenoid biosynthesis, diterpenoid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, zeatin biosynthesis, and brassinosteroid biosynthesis, indole alkaloid biosynthesis as well as plant hormone signal transduction were all enriched in the mutant ( Fig. 4 and Additional file 8: Table S5). In addition, MAPK signaling pathway was also enriched. These results suggested that plant hormone play major roles in regulation of the mutant phenotype. There were other pathways including photosynthesis, phenylpropanoid biosynthesis, flavonoid biosynthesis, isoflavonoid biosynthesis, and flavone and flavonol biosynthesis pathways were found to be enriched.
In leaf, brassinosteroid biosynthesis, indole alkaloid biosynthesis, plant hormone signal transduction and MAPK signaling pathways were enriched in accordance with those in stem ( Fig. 4 and Additional file 8: Table S5). Phenylpropanoid biosynthesis, flavonoid biosynthesis, isoflavonoid biosynthesis and plant-pathogen interaction were also in the top 20 enriched pathways. Different from that in stem, some glucide and glucolipid metabolic pathways were enriched in leaf.

Cell wall related genes
RNA-seq results showed that several CesA genes and all CSL genes were down-regulated in stem. The expression of most expansin genes was down-regulated in sdm2 stem compared with those in FH1 ( Fig. 5 and Additional file 3: Table S1). Interestingly, the expression of CesA, CSL and expansin genes was up-regulated in leaf ( Fig. 5 and Additional file 4: Table S2). These results were coincided with the phenotype of sdm2 leaf, dark green in color, thicker, smaller than FH1 leaf.
Some XTHs were significantly up-regulated and some were down-regulated in sdm2 stem ( Fig. 5 and Additional file 3: Table S1). In leaf, the expressions of all XTHs were up-regulated, and some XTHs were 16-fold higher than those in FH1 ( Fig. 5 and Additional file 4: Table S2). Some genes encoding endo-1,4-β-glucanases (EGases) in stem and leaf were down-regulated.
Expression changes of hormone biosynthesis and signaling genes DWF4 was designated as CYP90B1 and proposed to catalyze steroid 22α-hydroxylation in BR biosynthesis [64]. A gene annotated as cytochrome P450 90B1 showed 73% identity with Arabidopsis DWF4 gene, which play key roles in BR biosynthesis. In stem, the expression of this gene was detected in FH1, while not detected in sdm2. In leaf, the expression of this gene in sdm2 was also lower than that in FH1. Previous studies confirmed that cytochrome P450 monooxygenase CYP90A1 (CPD) acted as BR C-3 oxidase [36]. The expression level of CYP90A1 (CPD) in sdm2 stem decreased by 1.2-fold than that in FH1. While in leaf, its expression was slightly higher in sdm2 than in FH1. The cytochrome P450 family members, CYP90C1/D1, CYP85 (Dwarf) and CYP71A1, mediate rate-limiting reactions in BR biosynthesis [65][66][67]. Our results showed that the expression of CYP90C1/D1 gene decreased by 2.64-fold and 1.41-fold in stem and leaf respectively in sdm2. The Dwarf genes which encoding cytochrome P450 85A-like in stem Fig. 4 Bubble diagram of top 20 enriched KEGG pathways of DEGs in stem and leaf. X axis represents the Rich Ratio, which meaning the ratio of selected gene number annotated to a particular item to the total number of genes in this item in one species. The calculating formula is Rich Ratio = Term Candidate Gene Num/Term Gene Num. Y axis represents KEGG Pathway. The size of the bubbles indicates the number of genes annotated to a KEGG Pathway. And the color represents Q-value of enrichment. The deeper the color, the smaller the Q-value and cytochrome P450 85A in leaf were all up-regulated. CYP71A1 which catalyze the generation of deoxocastasterone (6-deoxoCS) and castasterone (CS) was downregulated both in stem and leaf in sdm2. A cytochrome P450 monooxygenase (CYP734A1) which degrades BRs was up-regulated in sdm2 leaf ( Table 2, Additional files 3 and 4: Table S1 and S2).
In both stem and leaf, the expression levels of BRI1 and BAK1 were down-regulated in sdm2. Noticeably, the expression levels of BES1/BZR1 gene were extremely low in stem and leaf in sdm2, 170-fold and 79-fold lower than in FH1, respectively. The expression levels of IWS1 in stem were up-regulated in sdm2, while in leaf, it was not differentially expressed. The expression of cyclin-D3-3 encoding gene was down-regulated in sdm2 stem, which was coincident with the expression trends of BR-related genes. While the expression levels of cyclin-D1 and cyclin-D5 increased in sdm2 leaf (Table 2, Additional files 3 and 4: Table S1 and S2).
Several GA biosynthesis and signal transduction genes were differentially expressed between sdm2 and FH1. In sdm2, most of GA20ox and 2-oxoglutaratedependent dioxygenase encoding genes, which promote GA biosynthesis, were up-regulated in stem, while most 2-oxoglutarate-dependent dioxygenase decreased in leaf. While two genes encoding GA3ox were down-regulated in stem and its differential expression was not detected in leaf. GA2ox convert excess GAs to inactive forms through 2β-hydroxylation [29]. In mutant stem and leaf, genes encoding GA2ox1 and GA2ox2 were decreased in expression levels while the expression of three genes encoding GA2ox8 were up-regulated in stem ( Table 2, Additional files 3 and 4: Table S1 and S2).
Two GID1 genes in stem showed opposite change trend and one GID1 gene was down-regulated in leaf. SCL3 (Scarecrow-like 3) seemed to attenuate DELLA protein and acted as a positive regulator in GA pathway [68]. In mutant stem, SCL3 gene was up-regulated while two SCL14 and one SCL21 genes were down-regulated. Interestingly, all genes encoding gibberellin-regulated protein (GRP) were reduced in sdm2 stem while upregulated in leaf, which were consistent with the change trends of cell wall related genes including CesA, CSL, expansin and XTH genes (Table 2, Additional files 3 and 4: Table S1 and S2).
Four AUX1 genes in stem and two in leaf were all significantly up-regulated, while the expressional levels of two genes encoding auxin efflux carrier (PIN) in leaf decreased. In addition, the expression change trends of MDR/ABCB genes in stem and leaf of sdm2 were irregular, most of which were downregulated and some were up-regulated. In sdm2 stem, two genes encoding AUX22 were down-regulated and almost all ARF members were up-regulated, while in leaf, the expressions of all Aux/IAA and ARF genes were increased. In addition, most SAUR genes in stem and almost all SAUR genes in leaf were up- Fig. 5 The expression pattern of cell wall related genes in stem and leaf of sdm2. Heatmaps represent the differential expression of CesA and CSL a, XTHs b, expansions c in stem and leaf between sdm2 and FH1. Scale bar is located at upside with log2 ratio value varying from green to red

Transcription factors involved in plant growth and development
Our results showed that several transcription factor genes involved in hormone signaling pathways and their downstream changed significantly. Gene encoding SHORT INTERNODES-like protein was downregulated both in stem and leaf of sdm2. Whereas, SHI-RELATED SEQUENCE 3-like (SRS) gene increased in stem and decreased in leaf. All WRKY30/ 46/70 genes were down-regulated in sdm2 stem and leaf. The expression levels of four PRE6 genes decreased in stem and three PRE6 genes were up-regulated in leaf ( Table 2, Additional files 3 and 4: Table S1 and S2).

Verification of DEGs using qRT-PCR
In order to validate the RNA-Seq data, the expression levels of DEGs were detected using Quantitative Realtime PCR (qRT-PCR). A total of 10 genes related to hormone signal transduction and cell wall organization were selected ( Fig. 6 and Additional file 9: Table S6). The relative expression levels (log 2 sdm2/FH1) of these genes estimated by qRT-PCR were generally consistent with those by RNA-seq. The overall correlation coefficient of a liner regression analysis in stem and leaf was 0.7737 and 0.8328, respectively (Fig. 6).   [4-7, 20-23, 40, 50, 69]. Microscopic assays indicated that the dwarf phenotypes of BR-deficient and BRinsensitive mutants were caused by reduced cell size [23,32]. Maize br2 mutant, deficient in the polar auxin transport, showed reduced stalk cell length and diameter [69]. In our study, sdm2 phenotype was resulted from reduction both in internode number and internode length. The difference might be due to changes in cell division and cell expansion during the stem growth process, which was supported by the identification of a number of cell wall biosynthesis and metabolism enzyme encoding genes as DEGs. Cell division and cell expansion need the biogenesis or remodeling of cell wall. The extensibility of cell wall controls the rate of cell expansion. Genes, including CesA and CSL, responsible for cell-wall biogenesis, and expansins and XTHs, involved in cell wall loosening, are all critical in plant growth and development. Our transcriptome data revealed several CesAs, almost all CSLs and expanin genes were downregulated in sdm2 stem. Some XTHs were downregulated in sdm2 stem. Besides, genes encoding EGases which functions in the modification of hemicellulose network were down-regulated in sdm2 stem. These results suggested that the cell wall biosynthesis and extensibility were reduced in sdm2 stem, which might result in reduced cell division and cell expansion and further cause short internode and small number of stem nodes. Studies revealed that many cell wall organization enzymes, including CesA, XTHs, EGases, and expansins were regulated by BRs [24,45,48,[70][71][72]. The essential roles of BRs in stem elongation have been confirmed by several mutants deficient in BR biosynthesis or perception [23,[30][31][32]40]. Our results showed that key BR biosynthetic genes CPD, DWF4 and CYP90C1/D1 and CYP71A1 were all down-regulated in sdm2. CYP734A1, which degrades BRs, was up-regulated in leaf of sdm2. The changing trends of above genes suggested the contents of BRs might be decreased in sdm2. BR response genes including BRI1, BAK1 and BES1/BZR1, were down-regulated, especially BES1/BZR1 in sdm2 stem. These results suggested that BRs could be a key factor in the formation of semidwarf phenotype of sdm2. In addition, BR signaling pathway also regulate the leaf and seed size and shape. In Arabidopsis cpd mutant, the decrease of cell size and cell number resulted in reduced leaf size [73]. The seeds of Arabidopsis det2 mutant were smaller and shorter because of decreased seed cavity, reduced endosperm volume and cell length and delayed embryo development [74]. Transgenic Arabidopsis and tobacco lines overexpressing DWF4 displayed increased plant height, longer petiole and leaf blade length, increased number of branches and siliques compared to the control [75]. In this study, sdm2 plants showed smaller leaf size, shorter petiole, and both decreased seed size and weight, which were all in accordance with phenotypic characteristics of BR-deficient mutants. We speculate that CPD, DWF4, BRI1 and BES1/BZR1 are key genes that lead to the phenotypic differences between sdm2 and FH1 line.
In this study, several GA biosynthetic genes showed irregular change trends and some seemed to promote the biosynthesis of GAs, while all GRP genes decreased in sdm2 stem which had the same expressional pattern with those of cell wall related genes. This result could be a reflection of lowered GA level in sdm2 stem, and the up-regulated expression of GA20ox and 2-oxoglutaratedependent dioxygenases might be resulted from the feedback regulation of GA signal pathway. The polar transport of auxin was essential for plant growth and development. In sdm2, the increased transcription of auxin influx carrier and down-regulation of auxin efflux carrier together with the irregular changes of auxin polar transporters, MDR/ABCB genes, suggested the disordered auxin transport from shoot apex and young leaves to stem. It might contribute to the formation of short internode in sdm2.
BRs could affect GA biosynthesis through positively regulating GA20ox expression [76]. DELLA protein could directly interact with BZR1 and inhibit its DNA binding ability, and the promotion of GAs on cell elongation require BZR1 or BRs [57]. Application of auxin could induce the expression of GA20ox and GA3ox while reduce the GA2ox transcript [77,78]. Auxin and BRs also could interact with each other and synergistically regulate plant development [79,80]. Auxin could directly induce DWARF4 expression and BR biosynthesis [81]. In return, BRs influenced auxin redistribution through stimulating the expression of PIN genes [82]. There are complicated cross talk among hormones during plant growth and development. Our results suggested that the significantly down-regulated BR biosynthetic and response genes, low GA response and disordered auxin transport acted synergistically and contributed to the defects in sdm2 stem growth.
As the common target genes of the core transcription module, PRE family factors connect external and endogenous signals with downstream cell elongation components together through an antagonistic cascade reaction. The expressions of PRE6 genes in sdm2 stem and leaf were consistent with those of cell wall enzymes, which further proved the direct promotion effects of PRE on cell elongation. Tracing back to the upstream hormonal signals, the extremely low expression of BR biosynthetic and signal transduction genes and low GA signaling in sdm2 stem all contributed to the down-regulation of PRE genes, and subsequently reduced cell elongation factors which finally caused the semi-dwarf phenotype. Based on these results, we proposed a potential model regulating peanut stem development (Fig. 7). BRs, GAs and auxin signaling pathways interact cooperatively to influence PRE expression, and regulate downstream cell wall-related genes to control cell elongation during stem development. However, further research is required to verify this hypothesis after cloning the mutant gene.

Conclusions
The transcriptome analysis of sdm2 identified a number of differentially expressed genes involved in hormone biosynthesis, signaling and cell wall synthetic and metabolic pathways. Especially several genes in BR pathway were significantly down-regulated in stem and leaf of sdm2 as compared to FH1. Genes in cell wall synthetic and metabolic pathway which related to cell elongation were generally down-regulated in sdm2 stem. These findings provide critical information for uncovering the molecular genetic control of peanut stem development.

Plant materials
The cultivated peanut cultivar Fenghua1 (FH1) was bred by and acquired from Yongshan Wan's laboratory in Shandong Agricultural University, Shandong Province of China. In 2013, the seeds of FH1 was irradiated with 60 Co γ-ray of 500 Gy. In M 2 generation, a line with phenotypic segregation of normal and semi-dwarf plant height was screened. After self-crossing for another two generations, the stable semi-dwarf peanut line was obtained and named semi-dwarf mutant 2 (sdm2).

Planting conditions and sampling methods
In 2018, sdm2 and its wild type FH1 were sowed in experimental farm of Shandong Academy of Agricultural Sciences. The ridging mode was adapted with ridge spacing of 90 cm, ridge width of 55 cm and ridge height of 12 cm. Peanut seeds were sowed for two rows on each ridge with a single seed per hole, and the plant spacing and row spacing were 18 cm and 35 cm, respectively. All field management followed the standard agricultural practices. For each line, about 100 plants were planted. It was sowed on May 1st and harvested on September 5th with the growing period of 129 days. The days after planting (DAP) was used to calculate seedling growth time. During the vegetative period, at least 5 plants from each line were selected to investigate the main stem height through measuring the length from base of the aboveground plant to the tip of the main stem for every 2 weeks. The 1~5 internodes and the first fully developed leaf were collected from 60 DAP plants of sdm2 and FH1 for RNA sequencing. At this time, the main stem height of sdm2 and FH1 was significantly different. The samples were frozen immediately in liquid nitrogen, and stored at − 80°C for RNA extraction. For both stem and leaf, three replicates were prepared.
The hydroponic culture experiment was conducted in illumination incubator in the laboratory. The temperature was 32°C and the light-dark cycle was 14 h of light and 10 h of darkness. Seeds of sdm2 and FH1 were planted in glass ware filled with Hogland culture solution with the preservative film for support. About 20 plants were planted for each line. At the 14 DAP, the number of internodes and the length of each internode and hypocotyl were measured. The main stem height was the summation of each internode length. Each line had five replicates.

RNA extraction, cDNA synthesis and high throughput sequencing
Total RNA was extracted from stem and leaf using Trizol Reagent (TaKaRa, Inc., Dalian, China) according to the manufacturer's instructions. RNA samples were treated with DNase I to remove genomic DNA contamination. RNA quality and quantity were analyzed using Agilent 2100 and NanoDrop. mRNA was enriched using magnetic beads with Oligo (dT) and cleaved into short fragments (~200 nt) in fragmentation buffer. The reverse transcription was conducted with random hexamer primer and then the second strand cDNA was synthesized. After end repair, the 5′ tails were phosphorylated, the 3′ tails were added with anadenine. Sequencing adaptors were ligated to the double-stranded DNA fragments. Then the fragments were amplified by PCR to construct cDNA library. The library was sequenced using BGISEQ-500 platform by Beijing Genomics Institute (BGI).

Bioinformatics analysis of RNA-Seq data
Raw reads were generated from each cDNA library. To obtain clean reads, the adaptor sequences, low-quality reads and reads containing more than 5% unknown bases were removed using SOAPunke and trimmomatic software. All clean reads were mapped to the reference genome of Arachis hypogaea cv. Tifrunner (https://www. peanutbase.org/data/public/Arachis_hypogaea/) using HISAT2 program [83]. The clean reads were aligned to reference gene by Bowtie2 (RNA-Seq by Expectation Maximization) [84]. The statistical analyses of randomness, degree of coverage, and sequencing saturation were also accomplished. The gene expression level was calculated with RSEM method [85] and normalized to FPKM (Fragments Per Kb per Million reads). The relative gene expression level between two samples was counted by log 2 ratio. Pearson's correlation coefficient between every two samples and principal component analysis (PCA) were performed. Differentially expressed genes (DEGs) were identified using DEGseq2 method and screened with the criteria of fold change≥2 and Q-value≤0.001 [86].
Gene Ontology (GO) annotation was carried out by Blast2GO program through comparing DEGs with GO terms in the GO database and GO functional classification was performed using WEGO software. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was conducted by mapping DEGs to KEGG database. The GO and KEGG functional enrichment were performed using hypergeometric test. The p-value formula in hypergeometric test can be acquired from https://en.wikipedia. org/wiki/Hypergeometric_distribution in detail. FDR (False Discovery Rate) correction of all p-values was conducted. The GO terms and KEGG pathways whose FDR ≤ 0.01 were defined as significant enriched. BLASTX (E < 0.00001) against NCBI Nr database was carried out. For transcription factor (TF) annotation, the ORF of DEGs were detected by getorf software, then these ORFs were aligned to TF protein domain (data from PlantTFDB) using hmmsearch program, and finally TFs were identified according to the TF family characteristics described in PlantTFDB (http://plntfdb.bio.uni-potsdam.de/v3.0/). The DEGs were also aligned to Plant Resistance Gene Database (PRGdb, http://prgdb.crg.eu/) by DIAMOND software to identified the resistance genes according to the query coverage and identity. The heatmap of cell wall related gene expression was charted using MeV software (https://sourceforge.net/projects/mev-tm4/) with the data of relative gene expression level between two samples counted by log 2 (sdm2/FH1).

qRT-PCR validation of RNA-Seq data
We used qRT-PCR method to verify the expressional levels of 10 selected genes. RNA samples were those used for high-throughput sequencing and the reverse transcriptions were performed using PrimeScript II 1st Strand cDNA Synthesis Kit (TaKaRa). The gene-specific primers were designed using PerlPrimer software and were listed in Additional file 10: Table S7. We performed qRT-PCR reaction on ABI7500 Real Time System (Applied Biosystems) using TB Green™ Premix Ex Taq™ II (TaKaRa). The parameters of thermal cycle were 94°C for 10 min, followed by 40 cycles of 94°C for 15 s and 60°C for 1 min in a 20 μl volume. Each reaction was performed three biological replications with actin gene as internal reference gene. The relative expressional level of each gene between sdm2 and FH1 was calculated by 2 -△△Ct method.
Additional file 1 Figure S1. The correlation heatmap of each sample of sdm2 and FH1. The X and Y axes represent each sample. The color