Response of peanut Arachis hypogaea roots to the presence of beneficial and pathogenic fungi by transcriptome analysis

Entomopathogenic fungus Metarhizium anisopliae obtain survival benefit meanwhile promote the nutrient absorption of root as an endophyte. However, little is known concerning molecular mechanisms in the process. We performed the transcriptome sequencing of A. hypogaea roots inoculated M. anisopliae and pathogenic Fusarium axysporum, respectively. There were 81323 unigenes from 132023 transcripts. Total 203 differentially expressed genes (DEGs) respond to the two fungi, including specific 76 and 34 DEGs distributed respectively in M. anisopliae and F. axysporum treatment. KEGG pathway enrichment for DEGs showed the two top2 were signal transductions of plant-pathogen interaction and plant hormone. By qRT-PCR, the mRNA level of 26 genes involved in plant-fungus interaction confirmed the reliability of the RNA-Seq data. The expression pattern of the key DEGs on jasmonic acid (JA) or salicylic acid (SA) signaling pathway presented regulating consistency with JA or SA concentration detected by HPLC-MS. Those significantly stronger down-regulated DEGs by M. anisopliae thanby F. axysporum linking to hypersensitive response and negative regulation of defense, and those specific up-regulated genes in M. anisopliae treatment may predict that the less immunity is conducive to symbiosis F. axysporum may trigger JA-mediated defense regulated by ERF branch of JA signaling pathway, whereas M. anisopliae does not.

Peanut Arachis Hypogaea is a globally important crop for food and oil extraction. It is often attacked by a range of insect species and diseases causing serious yield losses in cultivation. It is widely grown in the semi-arid tropics, wherein China contributes the highest share and India ranks second by 41.6% and 12.5% in world production, respectively. The soil-dwell insect such as root-gnawing white grubs Holotrichia spp. and root rot pathogen Fusarium oxysporum Schlecht. could critically damaged the root system and impact nut production [1][2][3] . Pesticides and biopesticides were required repeat applications for controlling the pests and diseases to protect yields. The entomopathogenic fungus Metarhizium anisopliae (Metchn.) Sorokin is widely used as a biocontrol agent to reduce crop damage by pests, and has shown high potential efficiency in control of soil-dwelling pests e.g. beet root maggot Tetanops myopaeformis, wireworms Agriotes spp and white grubs Holotrichia parallela [4][5][6] . While having a recognised role in insect control, M. anisopliae became of increasing interest because of its beneficial role on plant growth. It has been added to a list of fungi as a plant growth-promoter by rhizosphere competence and plant endophytes 7 . When M. anisopliae was applied to a cabbage experimental field at a rate of 10 13 spores per ha, the fungal density could remain at 10 5 propagules/g in the inner rhizosphere, while the amount was only of 10 3 propagules/g in nonrhizosphere soil after several months 8 . Tomato plants treated with M. anisopliae had significantly greater plant height, root length, and shoot and root dry weight than those of the untreated control, although the response depended on isolate and inoculation rate 9 . Growth promotion following inoculation with M. anisopliae was also found for switchgrass (Panicum virgatum), haricot beans (Phaseolus vulgaris), corn (Zea mays) and peanut (Arachis. hypogaea) [10][11][12]  Functional classification by GO, COG and KEGG. Gene Ontology (GO), an internationally standardized gene functional classification system, was used to classify the function of the predicted A. hypogaea unigenes. In total, 23299 unigenes with BLAST matches to known proteins were classified into three major functional ontologies using 1534 functional terms (Fig. 1a, Supplementary Table 1). As shown in Fig. 1a, the majority   Table 1). However, within each of the three major categories, a small number of genes were assigned to subcategories such as extracellular matrix part, channel regulator activity and guanyl-nucleotide exchange factor activity. In order to predict and classify possible functions, all unigenes were aligned to the Cluster of Orthologous Groups (COG) database in which orthologous gene products were classified. 10811 non-redundant unigenes were subdivided into 26 COG classifications (Fig. 1b, Table 2). The cluster related to general function prediction only (2108, 19.50%) was the largest group, followed by those for post-translational modification, protein turnover, chaperon (1404, 12.99%) and signal transduction mechanisms (948, 8.77%).
The Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway database is a knowledge base for the systematic analysis of gene functions in terms of networks of genes and molecules in cells. Based on KEGG analysis, 9841 unigenes were assigned to 267 pathways ( Table 2, Supplementary Table 2). The pathways involving the largest number of unique transcripts were signal transduction (962, 6.50%), which may be involved in fungi-plant interactions, followed by carbohydrate metabolism (950, 6.42%), whereas signaling molecules and interaction (8) was the smallest group (Fig. 1c, Supplementary Table 2).
Global differential expression genes. To reveal the molecular mechanism for plant and fungi interactions, the differential expression genes (DEGs) induced by inoculation with M. anisopliae and F. axysporum were analyzed. The resulting Pearson's correlation coefficient (R 2 ) between samples were quite high in AM vs AC (R 2 = 0.843), AF vs AC (R 2 = 0.843) and AM vs AF (R 2 = 0.861) ( Supplementary Fig. 2). The expression level of each assembled transcript sequence in different samples was measured through RPKM (Reads per kilo-base per million reads) values, and DEGs (q-value < 0.005 and log2 (fold change) > 1) were defined as genes that were significantly enriched or depleted in one sample relative to the other.
Based on the criteria above, of the 81323 unigenes, 164, 122 and 20 were detected as significantly different in the comparison between AM vs AC, AF vs AC and AM vs AF, respectively. Among these DEGs, 16 were up-and 148 down-regulated, respectively in AM vs AC. For AF vs AC, 28 were up-and 94 down-regulated, respectively. For AM vs AF, 13 up-and 7 down-regulated, respectively ( Supplementary Fig. 3). As shown in the Venn diagram, a total of 203 DEGs were shared by AM, AF and AC. Individual inoculation with pathogenic fungus F. axysporum specifically regulated 34 genes. 76 DEGs with special expression pattern in AM vs AC probably functioned combinative symbiotic process (Fig. 2).
Functional distribution and connection of differentially expressed genes. All the DEGs in inoculating treatment with M. anisopliae (AM) and F. axysporum (AF), were hierarchically clustered into eight sub-clusters based on k-means and identified. For the two fungal treatments compared to the control (AC), the genes with similar down-regulation pattern were assigned to subcluster 3, 5, 6 and 7, and the genes showing similar up-regulation were assigned to subcluster 2, 4 and 8. The genes in subcluster 2 showed significantly higher up-regulation in AF compared to AM, while the genes in subcluster 4 showed significantly higher up-regulation in AM compared to AF. In particular, the genes in subcluster 1 were up-regulated in AM but down-regulated in AF (Fig. 3).
Further, the functional distribution of the up-and down-regulated DEGs were compared in different treatments. Based on GO annotation, only 16 genes in total DEGs were significantly enriched on the two GO stems (Corrected p-value < 0.05) ( Table 4). They were all appeared in the M. anisopliae treatment, with nucleic acid binding transcription factor (GO:0001071) and sequence-specific DNA binding transcription factor (GO:0003700) activity. The directed acyclic graph by topGO showed that sequence-specific DNA binding transcription factor activity (GO:0003700) was subject to nucleic acid binding transcription factor activity (GO:0001071), and then both were classified into the main stem of molecular function (GO:0003674). The 16 DEGs of peanut response to M. anisopliae included a variety of transcription factor, such as WRKY, MYC, TGA,  Table 3. Genes involve in plant-fungus interaction.
ethylene-responsive transcription factors, dehydration-responsive element-binding proteins and nitrate transporters (Table 4). Biological behaviors of DEGs are complicated, being related to many inner-cell metabolic pathways. KEGG pathway enrichment analysis for DEGs revealed the combinative pathways that multidirectional connecting to defense-linked regulation. As showed in Table 5, the interesting results were that signal transductions of plant-pathogen interaction and plant hormone were the top two enriched pathways for the down-regulated DEGs in both AF vs AC and AM vs AC (Corrected p-value < 0.05). They went through MYC or ERF branch on the jasmonic acid (JA) signaling pathway and defense-linked RPM1, CDPK and WRKY points. Individually, glycine, serine and threonine metabolism and phenylalanine metabolism were distinctly enriched for the up-regulate DEGs in AF vs AC (Corrected p-value < 0.05), while nitrogen metabolism, lysosome degradation and antigen processing and presentation pathways were distinctly enriched for the up-regulate DEGs in AM vs AC (Corrected p-value < 0.05). Particularly, nitrogen metabolism, amino sugar and nucleotide sugar metabolism were enriched in AM vs AF DEGs (Corrected p-value < 0.05). The regulating change converged from multi-directions was able to correct some metabolic pathways.
In spite of the intricate metabolic pathways, a limited signal crossing or gene interacting still could be found and speculated. KEGG pathway showed that phenylalanine metabolic pathway was linked to salicylic acid (SA) biosynthesis. Also, abscisic acid (ABA) biosynthesis was showed linking to carotenoid biosynethesis pathway. L-phenylalanine was known as a substrate for the synthesis of SA, besides that, it was an influence factor on carotenoid aggregation. Although nitrogen metabolic pathways were enriched in both AM vs AC and AF vs AC, nitrogen-transporting protein appeared significantly high, along with nitrate reductase in AM vs AC. In contrast, some immune-related metabolic pathways such as lysosome degradation and antigen processing and presentation were enriched in AF vs AC ( Table 5).
Detection of salicylic acid (SA) and jasmonic acid (JA) in roots of peanuts under different treatments. HPLC-MS was used to detect the concentration of JA and SA in roots of peanuts, which were separately treated by distilled water (AC), F. axysporum (AF) and M. anisopliae (AM) (Fig. 5). The results showed that  there was no significant difference between AM and AC in the concentration of JA, which was around 320 ppm, while the concentration of JA in AF was significantly higher than that in AC or AM, which getting 750 ppm, The concentration of SA seemed similar among AC, AF and AM, which ranged 361-486 ppm. They were not significantly different between AC and AF and between AF and AM. However, the concentration of SA in AM was significantly lower than that in AC. The statistically significant difference of hormone concentration was considered on an error probability of p < 0.05 by one-way ANOVA (SPSS version 16.0, SPSS) using Duncan's multiple-range test.
Comparing to the DEGs in JA and SA signaling pathways that verified by qRT-PCR as above, we found that the concentration of JA was highest in AF but relatively low in AC and AM, which consistent with the highest expression level of defensin (c38108) by qRT-PCR as well the highest mRNA level of defensin by transcriptome analysis in AC. It is known that defensin is a marker gene in ERF branch of JA signaling pathway. Another coincidence is that mRNA level of MYC transcript factor, a key gene on the JA signaling pathway, is also lower in AC and AM. On the other hand, in SA signaling pathway, mRNA level of those DEGs including transcription factors WRKY (c50433, c48617), MYB (c44072) and TGA (c46800) and pathogenesis-related protein (c40732) are lower in AF and AM than in AC, that consistent with the SA concentration in the above detection. The results provided one more support for RNA-Seq data reliability.

Discussion
Based on the regulated direction, function distribution and KEGG pathway enrichment of the DEGs in AM vs AC and AF vs AC, a comprehensive integration may resolve the early different response to beneficial M. anisopliae and pathogenic F. axysporum in peanut roots, and provide clues to mechanisms underpinning specific symbiosis by M. anisopliae. Similar response of peanut root to inoculation of beneficial and pathogenic fungi. Of the 203 differentially expressed genes (DEGs) identified in AM vs AC and AF vs AC comparisons, 192 (c.94.6%), appeared in the same direction of up-or down-regulation (Fig. 3). It suggested that there were similar pattern in plant cellular response and interaction with pre-invading fungi, whether they are beneficial or pathogenic to plant. The pattern may resemble or equivalent to the first tier of plant immune system. Plant evolved pattern recognition receptors (PRRs) can be activated by recognition of evolutionarily conserved pathogen-associated molecular patterns (PAMPs) which then trigger the mitogen-activated protein kinase (MAPK) cascades. In turn, this activates plant hormone signaling to integrate various aspects of the multi-layered plant defense response (called pattern-triggered immunity or PTI) 16,17 . The optional signaling may guide different types of immune reaction 18 . Furthermore, it is also acknowledged that fungi can produce effectors to suppress the basal PTI by suppressing immune signaling, inhibiting MAPK cascade or blocking resistant protein expression 19,20 . These co-evoluted effectors are deployed to modify host cell processes or to associate with the host so that reaching biotrophic lifestyle for part, if not all, of their lifecycle 21 . In turn, plant disease resistance proteins (R proteins) could recognize pathogen effectors to induce effector-triggered immunity (ETI) which usually cause hypersensitive response (HR) 22,23 . Colonization of the roots by either M. anisopliae or F. oxysporum may mean that they are recognized by the same PRRs located on the surface of root cells, which consequently leads to similar cellular responses as found in this study. In the transcriptome of this experiment, the part of genes up-regulated may be activated to identify and defend against invading pathogens. We found the common up-regulated genes from subclusters 2, 4 and 8 in Fig. 3, linking to several biological processes including chitinase, pro-hevein, peroxidase, proteinase, nitrate reductase and lectin. These proteins and enzymes are directly or indirectly involved in the regulation of plant defense. For example, vacuolar processing enzyme (VPE), (c49130_g1) which was up-regulated in AM vs AC and AF vs AC, may be a cysteine protease responsible for caspase -1 activity and promote programmed cell death (PCD) by disrupting the vacuole in pathogenesis like in tobacco 24 . It may positively regulate HR cell death mediated by NB-LRR recognition of the invading pathogen like metacaspase-1 in a variety of plants 25 . The up-regulated 1-amino-cyclopropane-1-carboxylate oxidase (ACO) may be involved in ethylene synthesis, reasoning from 1-amino-cyclopropane-1-carboxylic acid (ACC) is ethylene precursor 26 . Lectins are a group of sugar-binding  proteins which could regulate immune response to guide plant defense 27 . Likewise, chitinases often contain a hevein domain which related to defense by carbohydrate degradation 28,29 . Besides, one actin-like pathogenesis related protein and one defensin also found in up-regulated DEGs in AM vs AC and AF vs AC.
Another part of genes, down-regulated in both AM vs AC and AF vs AC, may be suppressed by the two fungi to facilitate subsequent invasion. We found commonly recorded down-regulated genes from subcluster 3, 5, 6 and 7 in Fig. 3, that referred to a variety of physiological processes, such as plant-pathogen interaction, plant hormone signal transduction, amino acids metabolism, starch and sucrose metabolism. The similar response to both fungi may be required to directly deal with fungal invasion or to trigger subsequent immune reaction in peanut cells. For example, we identified nine TIR-NB-LRR or CC-NB-LRR resistance proteins, as well as one LRR receptor-like serine/threonine-protein kinase, one chitin elicitor receptor kinase and three receptor-like protein kinase. They may play important role in recognizing specific effectors. In regards to SA, a plant hormone closely related to plant disease resistance, we found a down-regulated heat shock protein HSP70 that was thought to process a defense-promoting function by promoting SA pathway 30 . Some enzymes related to SA synthesis such as one EDS1/PAD4, one tryptophan synthase and one arogenate dehydratase (choloroplastic-like), and some transcription factors increasing SA accumulate by positive regulation of ICS gene encoded in isochorismate synthase pathway, such as MYBs, WRKYs and WIPK, also were down-regulated. On the jasmonic acid signaling pathway, we found two pathogenesis-related marker genes PDF1.2 and VSP2, their down-regulation would be for MYC and ERF branches for JA signaling suppression 17 . The MYC branch and the ERF branch are two major branches of the JA signaling pathway 31 . There were four down-regulated MYC and ERF transcription factors in AM vs AC and AF vs AC, they may be repressed by jasmonate ZIM-domain (JAZ) transcriptional repressors. The COI1-JAZ co-receptor complex leads to ubiquitination and proteasome-dependent degradation of JAZ repressors and release of MYC or ERF transcription factors from transcriptional repression [32][33][34] . Ubiquitination regulates pattern-recognition receptor signaling that mediates immune responses 35 . There are four E3 ubiquitin protein ligase down-regulated in AMvsAC and AFvsAC. On the other hand, we found down-regulated (+)-abscisic acid 8′-hydroxylase involved in the ABA catabolic pathway. ABA promotes MYC branch by protein phosphatase 2C (PP2C) and antagonizes ERF branch 31,36 . In addition, A chitin-inducible gibberellin-responsive protein was found down-regulated. It could regulate the growth-repessing DELLA proteins which acts positively on JA signaling by sequestering JAZ repressor 37,38 . DELLA proteins maybe play an important role in hormonal cross talk to repress defense responses in symbiosis 39 .
As described above, down-regulated DEGs may be suppressed to weaken immune and inhibit development. These multi-pathway process may be very complex, but they do provide us with some important joints resulting transcriptomic and metabolomic data subsets.

Specific response of peanut roots to the presence of beneficial fungus compared to pathogenic fungus. The genes, up-regulated expression in M. anisopliae treatment but down-regulated expression in
F. axysporum treatment, should be particularly concerned (Fig. 3, subcluster 1). Ferric Redictase Defective 3 (FRD3), a MATE family member, has been shown to be an efflux transporter of the efficient iron chelator citrate in A. thaliana. It mediates-citrate release in the apoplastic space, and represents an important process by which efficient iron nutrition is achieved between adjacent tissues lacking symplastic connections, maintaining iron homeostasis throughout plant development 40,41 . Iron was identified as an essential micronutrient for the legumerhizobium symbiosis 42 . Another special DEG c50341_g1 up-regulate in AM vs AC but down-regulate in AF vs AC, has acetyltransferase activity and a NB-ARC domain which makes us surmise it as a receptor to symbiotic effectors. And more, some lipid transfer proteins should be required for the successful symbiotic association between a microbe and its host, for example, lipid transfer MtN5 for the symbiosis of Sinorhizobium meliloti and Medicago truncatula 43 . The role of high affinity nitrate transporter which was significantly up-regulated in AM vs AC may be closely contacted with the evidence of M. anisopliae providing more nitrogen sources for the plant 6,44 . Mutual benefits were an important basis for symbiosis. The other two enzymes, Feruloyl CoA ortho-hydroxylase 1 and 3′-N-debenzoyl-2′-deoxy taxol N-benzoyltransferase, have not been reported as having a role in plant-fungal symbiosis. But there is a possibility that they have an important role on synthesis of scopoletin and taxol, respectively 45,46 . Scopoletin has been shown to have distinctly antifungal synergistic effects in Melia azedarach L. 47 . Taxol, usually extracted from plants and used as antitumor, was identified antimitotic activity produced by endophytic fungi 48 .
Furthermore, in the above 3.1 description, those DEGs in the same direction of up-or down-regulation presented different regulatory levels in M. anisopliae treatment and in F. axysporum treatment. After a detailed analysis on their functions and metabolic pathways, we were surprised to find that almost all of significantly stronger down-regulated DEGs in AM than in AF were linked to those transcription factors involved in hypersensitive response (HR) and in negatively regulating defense (Supplemental Table 1). HR represents a strong incompatibility and a rapid response to prevent the spread of microbial infection by PCD in plant. In this experiment, the HR factors CaM, CDPK, PBS, RPM1 and RPS were down-regulated, which may serve to weaken the strength of the defense response in the plants when interacting with the fungi. This phenomenon was also found in other pathogenic fungi infecting to plants. Barley powdery mildew Blumeria graminis and the leaf spot disease Stemphylium lycopersici can suppress HR-based PCD [49][50][51] . Phytophthora infestans and Mycosphaerella pinodes can suppress HR-based defense when they infect host plants 52,53 . We also found that all of those down-regulated resistance proteins and related kinase in AM vs AC and AF vs AC, including nine TIR/CC-NB-LRR, one LRR receptor-like serine/threonine-protein kinase, one chitin elicitor receptor kinase and three receptor-like protein kinase, were expressed significantly lower in AM vs AC than in AF vs AC. Integrated analysis of the transcription factors down-regulated found that most of them involved in negatively regulating defense by their respective paths, including JA-JAZ-ERF/MYC, FLS-MEKK-WRKY and GA-GID2-DELLA, or by cross regulation of JA path with SA, ET or ABA. Conversely, almost all of significantly stronger up-regulated DEGs in AF than in AM were assigned to enzymes and proteins involved in positively regulating defense. They include PR, peroxidase, chitinase, lectin and indole glucosinolate (Supplemental Table 1). Therefore, these results suggest that, even though both the two fungi encountered a plant resistance defense, the path and strength of the response were different. For M. anisopliae, the defense was weaker and came from indirect multipath of hormone regulation. The indirect multipath regulation may be accompanied by other responses as description above. For F. axysporum, the defense was stronger and more direct. It may be because F. axysporum triggered JA-mediated defenses against necrotrophic fungal pathogens, which regulated by ERF branch of JA signaling pathway, whereas the non-pathogenic fungus M. anisopliae did not 54 .
As it was well known, plant defense implemented by a refined immune system in which various hormones interact forming complex network. In spite of only 4 d interaction between peanut and the two fungi in the experiment, the DEGs distribution can still expose some clues of hormones cross-regulation. The transcriptional factors MYC and ERF were branched to two major lines of the JA signaling pathway 31 . The MYC controlled the branch including the downstream marker gene vegetative storage protein 2 (VSP2), while the ERF as an ethylene response factor regulated the other branch including the marker gene plant defensing (PDF) 55,56 . The ERF branch of the JA pathway is associated with enhanced resistance to necrotrophic pathogens. SA activate PRs (pathogen related genes) by the TGA transcription factor 57,58 . Heat shock protein 70 (HSP70) is thought to possess a defense-promoting function by promoting the SA pathway 30 . There is also cross talk between SA and JA pathway. The growth-repressing DELLA proteins are regulated by gibberellins which acts positively on JA signaling by sequestering JAZ repressor 37,38 . Ethylene acts synergistically on the expression of ERF branch of JA pathway, whereas it antagonizes the MYC branch 59 . However, ABA promotes MYC branch by PP2C and antagonizes ERF branch 31,36 . WRKY transcription factors are essential for plant defense response 60 . LRR receptor-like serine/threonine-protein kinase FLS2 activate WRKY25 and WRKY33 to repress defense-related gene induction by MAPK cascades signal transduction 61 . Recently, cross regulation of plant hormone was reported more interpretation 31 . In this study, it can be inferred that cross regulation of peanut hormones guided the down-regulated expression of numerous transcription factors in both M. anisopliae and F. axysporum treatments. Moreover, those transcription factors, significantly stronger down-regulated in AM than in AF, directly or indirectly resulted in changing activity of the enzymes and proteins related to defense by negative regulation. Oppositely, the enzymes and proteins related to defense such as PDF (plant defensin), PRs (pathogen related proteins), peroxidase, chitinase and lectin, etc. were significantly higher in AF than in AM. This indicates that in the early stages of infection, both beneficial M. anisopliae and pathogenic F. axysporum suppressed some important transcription factors by various hormones pathway and inspired plant resistance. The difference is that F. axysporum caused stronger immunity than M. anisopliae done. The predictable result is that the less immunity in M. anisopliae treatment could be conducive to establish a symbiotic relationship.

Supposing the pathway factors of M. anisopliae symbiosis in plant. Beneficial interaction could
help M. anisopliae colonization and symbiosis in plant. Based on the results and analysis of this study, we speculate that M. anisopliae symbiosis in plant might be divided into two parts. One part might contain weakening HR and regulating defense response, nutrient supply and some specific induction, which related to those significantly stronger down-regulated transcription factors in AM than in AF. Another part might be similar to known common symbiotic way via a central pathway consisted of symbiosis receptor-like protein kinase (SymRK), Ion channel (CASTOR/POLLUT) and calcium-and calmodulin-dependent protein kinase (CCaMK), as if azotobacter in legume and arbuscular mycorrhizal fungi in various plant. The latter part had not yet been induced and differently expressed in the early interactions of plant and fungi in this experiment. Another possibility is that M. anisopliae may have a different pathway to establish successful symbiosis.

Materials and Methods
Sample collection and preparation of peanut Arachis hypogaea L. roots. To investigate root response to pathogen invasion of beneficial and harmful fungi we established three treatments: A. hypogaea inoculated with M. anisopliae (AM), inoculated with F. axysporum (AF) and an un-inoculated control (AC). The treatments were established in aseptic condition to minimize cross contamination from other microorganisms associated with the A. hypogaea rhizosphere, including harmful and beneficial to roots. Furthermore, it is difficult to determine the progression of infection by infection time because M. anisopliae colonizes the rhizosphere for an indeterminate period before becoming endophytic and enter the root system. The testing peanut A. hypogaea species was the Chinese cultivar Luhua-11. The soaking seeds were sown in pots containing sterilized vermiculite, one seed in one pot. The plants grow in a controlled climate chamber at 25 °C and L:D 14:10. After seven days, when seedlings were at the two lateral branches stage, 20 mL conidial suspensions of M. anisopliae (1 × 10 7 spore/ mL) or of F. axysporum (1 × 10 7 spore/mL) were carefully drip irrigated to each root, respectively. The control plants received 100 mL of sterile water. Four days following treatment, the samples were collected by carefully pouring pots, sweeping aside vermiculite, removing the plants and washing the roots in sterile water. The roots were cut and immediately immersed in liquid nitrogen, then stored at −70 °C. To increase the accuracy of detection, three subsamples of peanut roots were bulked into one sample for each RNA seq sample.
RNA extraction, library construction and sequencing. Total RNA was isolated from frozen medullar tissue by using the RNA plant mini kit with column DNase digestion (Qiagen, Hilden, Germany) following the manufacturer's instructions. RNA degradation and contamination was detected on 1% agarose gels. RNA concentration was then measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). Additionally, RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
Scientific RepoRts | 7: 964 | DOI:10.1038/s41598-017-01029-3 A total of 3 μ g RNA per sample was used as input material for the RNA preparations. Finally, three samples with RNA integrity number (RIN) values above 8 were used for construction of the libraries. Sequencing libraries were generated using NEBNext Ultra ™ RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer. First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH − ). Subsequently, second strand cDNA synthesis was performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferential 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C. This was followed by PCR performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primers, respectively. Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and 125 bp paired-end reads were generated.
Sequence reads mapping, assembly and annotation. Raw data (raw reads) of Fastq format were firstly processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and the sequence duplication level of the clean data were calculated. All the downstream analyses were based on high quality clean data.
The left files (read1 files) from all samples were pooled into one big left.fq file, and right files (read2 files) into one big right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity with min_kmer_cov set to a default value of 2 and all other parameters set to default 62 .
Gene function was annotated based on the following seven databases: Nr (NCBI non-redundant Protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database) and GO (Gene Ontology). Data for each sequenced library was analysed using BLAST with a cutoff E-value of 10 −5 .
Differential expression analysis. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edge R program package through one scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq (2010) R package. The p-value was adjusted using q-value 63 . q-value < 0.005&|log2 (fold change)| > 1 was set as the threshold for significantly differential expression.
GO and KEGG enrichment analysis of differentially expressed transcripts. Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GO seq R packages based on the Wallenius non-central hyper-geometric distribution, which can adjust for gene length bias in DEGs 64 .
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 (http://www. genome.jp/kegg/) 65 . We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways 66 . Confirmation of the expression profiles by qRT-PCR. Differentially expressed genes identified by the above described method were validated using quantitative real-time PCR (qPCR). The real-time PCR was performed with the SYBR Premix ExTaqTM (TaKaRa, Dalian, China) on the ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The beta-tubulin gene was used as a reference control. The reaction was performed using the following conditions: denaturation at 95 ration 60 s, followed by 40 cycles of amplification (95 ol. The15s, 60 60 llo 60 s). Each plate was repeated three times in independent runs for all reference and selected genes. Gene expression was evaluated by the 2 −ΔΔCt method 66 .
Statistical analysis. For each sample, three technical replicates of the qRT-PCR assay were used. Results were expressed as means ± standard error (SE) of the number of experiments. The statistically significant difference of gene expression was considered only on an error probability of p < 0.05 by one-way ANOVA (SPSS version 16.0, SPSS) using Duncan's multiple-range test.
Gradient elution was performed withmobile phase A (0.1% formic acid in water) and B (acetonitrile) at 0.4 mL/min flow rate and at 40 °C. The initial composition (10% B) was maintained for 1 min, increased from 10% to100% B for 10 min, and returned to initial conditions over 23 min. A 1 min equilibration followed, yielding a total run time of 35 min.
The Q Exactive mass spectrometer was equipped with heated electrospray ionization source (HESI-II) and operated in the positive ionization mode. We optimized the following parameters: spray voltage 3700 V, capillary temperature 320 °C, heater temperature 425 425sheath gas and auxiliary gas flow were optimized at 30 and 10 arbitrary units respectively.
Preparation of LC-MS standard solutions and sample. Five concentrations (10 mg/L, 50 mg/L, 100 mg/L, 500 mg/L, 1000 mg/L) of SA and six concentrations (10 mg/L, 50 mg/L, 100 mg/L, 500 mg/L, 1000 mg/L, 5000 mg/L) of JA solutions were seperately used for the establishment of calibration curves. The stock solutions were stored at 4 °C. Calibration curves gave the respective equations y = 1699580*X, R² = 0.9956 and y = 672873*X, R² = 0.9965. 1 g of each peanut root samples were taken to grind the homogenate, and then add 1.5 mL 1:1 (v:v) mixture of ethanol and methanol to the homogenate. Ultrasonic for 2 hours followed by 13000 g centrifugation. The supernatant was stored at 4 °C for LC-MS detection. The sample solution was filtered through a 0.22 μm filter before LC-MS injection.
Results were expressed as means ± standard error (SE) of the number of experiments. The statistically significant difference of hormone concentration was considered only on an error probability of p < 0.05 by one-way ANOVA (SPSS version 16.0, SPSS) using Duncan's multiple-range test.