Skip to main content

Advertisement

We’d like to understand how you use our websites in order to improve them. Register your interest.

Integrated small RNA and mRNA expression profiles reveal miRNAs and their target genes in response to Aspergillus flavus growth in peanut seeds

Abstract

Background

MicroRNAs are important gene expression regulators in plants immune system. Aspergillus flavus is the most common causal agents of aflatoxin contamination in peanuts, but information on the function of miRNA in peanut-A. flavus interaction is lacking. In this study, the resistant cultivar (GT-C20) and susceptible cultivar (Tifrunner) were used to investigate regulatory roles of miRNAs in response to A. flavus growth.

Results

A total of 30 miRNAs, 447 genes and 21 potential miRNA/mRNA pairs were differentially expressed significantly when treated with A. flavus. A total of 62 miRNAs, 451 genes and 44 potential miRNA/mRNA pairs exhibited differential expression profiles between two peanut varieties. Gene Ontology (GO) analysis showed that metabolic-process related GO terms were enriched. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses further supported the GO results, in which many enriched pathways were related with biosynthesis and metabolism, such as biosynthesis of secondary metabolites and metabolic pathways. Correlation analysis of small RNA, transcriptome and degradome indicated that miR156/SPL pairs might regulate the accumulation of flavonoids in resistant and susceptible genotypes. The miR482/2118 family might regulate NBS-LRR gene which had the higher expression level in resistant genotype. These results provided useful information for further understanding the roles of miR156/157/SPL and miR482/2118/NBS-LRR pairs.

Conclusions

Integration analysis of the transcriptome, miRNAome and degradome of resistant and susceptible peanut varieties were performed in this study. The knowledge gained will help to understand the roles of miRNAs of peanut in response to A. flavus.

Background

Peanut (Arachis hypogaea. L), or groundnut, is one of the most important oil crops in the world. It is wildly planted in Asia, Africa, and North America. China, India and the United States are the world’s major producers of peanut. The world’s total volume of peanut production is about 29 million metric tons per year (http://www.peanutsusa.com), contributing to an estimated value of about $35 billion [1]. However, peanut production is vulnerable to the threat of aflatoxin contamination, which is caused by the infection of Aspergillus flavus and Aspergillus parasiticus [2]. The aflatoxin was first revealed in 1960, when the aflatoxin-contaminated peanut caused the death of a large number of turkeys in UK. Aflatoxins, can cause DNA damage and chromosome instability in host cells, are highly toxic and one of the most carcinogenic substances [3,4,5]. Peanut aflatoxin contamination is serious problem in China and limits the exportation of peanuts. Contamination can occur in the field, during harvest, or in storage and processing, and therefore, aflatoxin contamination is hard to prevent and control.

Although it is a widespread problem, the underlying molecular determinants and mechanisms of peanut-A. flavus interaction have remained elusive. The analysis of the molecular mechanism of peanut-fungus interactions, and mining peanut resistance information is the basis to develop resistance peanut varieties. To better understand the mechanisms, cDNA libraries derived from developing peanut seeds from a resistant cultivar (GT-C20) and a susceptible cultivar (Tifrunner) were constructed and a total of 21,777 EST sequences were generated [6]. Soon afterwards, a peanut oligonucleotide microarray chip was used to identify resistance genes to A. flavus and A. parasiticus mixed infection, and found that 62 genes in resistant cultivar were up-regulated and the expression of 22 putative Aspergillus-resistance genes were higher in the resistant cultivar in comparison to the susceptible cultivar [1]. In order to identify proteins related to the resistance to aflatoxin contamination, two-dimensional electrophoresis was employed with mass spectrometry as well as real time RT-PCR method and 12 potential different expressed proteins were identified [7]. Another group also used the differential proteomic approach and identified an array of proteins responding to A. flavus infection [8]. Recently, Wang et al. employed Next Generation Sequencing to survey the gene expression profiling of resistant and susceptible peanut genotypes in response to A. flavus growth [9, 10]. These studies provided valuable information for understanding the mechanism of peanut in response to aflatoxin production. However, due to the lack of peanut genomic sequences, it is hard for these studies to provide a comprehensive interpretation of the transcriptomic changes of peanut in response to A. flavus. Recently, the release of whole genome sequences of the diploid wild peanut species, A. duranensis and A. ipaensis [11, 12] and cultivars, Tifrunner and Shitouqi (https://peanutbase.org/, http://peanutgr.fafu.edu.cn/) provide new opportunity to understand the peanut resistance mechanism in response to A. flavus.

MicroRNAs (miRNAs) are a class of non-coding RNAs regulating gene expression and play roles in plant development and response to stresses [13,14,15]. Up to now, more than 5000 miRNAs have been identified from near 70 plant species (http://www.mirbase.org/). The functions of plant miRNAs in biotic stress responses have been extensively explored. For example, miR398 regulated plant productivity under oxidative stress conditions by targeting Cu/Zn superoxide dismutase [16, 17]. In tomato, the accumulation of miR169 was induced by drought stress. The expression of its target genes, nuclear factor Y subunit genes and multidrug resistance-associated protein gene, were significantly down-regulated by drought stress. Over-expression of miR169 enhanced drought tolerance of tomato [18, 19]. Many transcription factors, such as SPL [20, 21], MYB [22], NAC [23,24,25], ARF [26], HD-Zip [27] and AP2 [28, 29] were the targets of miRNAs, suggesting miRNAs located in the core position of the network in gene regulation.

Recently, the roles of miRNA in plant disease resistance were gradually discovered, and miRNAs were considered to play key roles in plant natural immune system [30, 31]. In Arabidopsis, it was proven that miR472/RDR6 silencing pathway modulated PAMP and effector triggered immunity (ETI) through the post-transcriptional control [32]. In rice, multiple miRNAs are involved in resistance against the blast fungus Magnaporthe oryzae. A group of rice miRNAs were differentially expressed upon M. oryzae infection. Over-expression of miR160a and miR398b enhanced resistance to M. oryzae in the transgenic rice [33]. Multiple miRNA target genes were related to disease resistance. For example, the expression and the variation of miR482 and miR2118 reflected the shift in the balance of plants NBS-LRR defense system [34]. In our previous studies, we identified many known and novel miRNAs in peanut using high-throughput sequencing method [35]. We found that the expression of some known and novel miRNAs was induced or inhibited upon Ralstonia solanacearum infection in cultivated and wild peanuts. Degradome sequencing results possibly indicated that some defense response genes were degraded by miRNA [36]. These data suggested that miRNA might play important role in peanut disease resistance.

In this study, we aimed to identify potential A. flavus responsive miRNAs and their targets in peanut, and study their roles in peanut-A. flavus interaction. For this purpose, we employed high-throughput sequencing technology to simultaneously analyze the miRNA and mRNA expression profiles and their regulation in response to A. flavus growth using a resistant and a susceptible peanut variety. Our results provided valuable information for understanding the molecular mechanism of peanut resistance to A. flavus.

Results

Different resistance of peanut varieties GT-C20 and Tifrunner to A. flavus

The resistance of peanut varieties GT-C20 and Tifrunner to A. flavus was validated. As a result, in two days after inoculation (DAI), hyphae were only observed in the seeds of Tifrunner, but not in GT-C20 (Fig. 1). From three DAI, the hyphae were very clear in both GT-C20 and Tifrunner. More hyphae were observed in the seeds of Tifrunner than that in GT-C20 (Fig. 1). Statistical analysis showed that GT-C20 was high resistant (Infection index of 4.7) and Tifrunner was high susceptible (Infection index of 95.3) to A. flavus growth according to the standard reported [37]. Samples of two DAI were collected for miRNAome and transcriptome analysis.

Fig. 1
figure1

Phenotype of peanut variety GT-C20 and Tifrunner after infected with Aspergillus flavus

High-throughput sequencing and identification of miRNAs

To study the role of peanut miRNAs in response to A. flavus growth, eight small RNA libraries were constructed and sequenced using Illumina HiSeq2000 platform. These small RNA libraries included two A. flavus-treated (GT1 and GT2) and two control (GC1 and GC2) libraries of resistant cultivar GT-C20, and two A. flavus-inoculated (TT1 and TT2) and two control (TC1 and TC2) libraries of susceptible cultivar Tifrunner. After sequencing, a total of 111,839,986 raw reads and 111,324,210 clean reads were obtained from eight libraries (Additional file: Table S1). The normalized clean reads were used for analysis of small RNA distribution. Most of the clean reads were 20–24 nucleotides (nt) in length, which accounted for 96.22, 96.67, 94.06 and 95.81% of all clean reads in TC, TT, GC and GT libraries, respectively (Fig. 2). The 24-nt small RNAs was the most abundant group, which was followed by the 21-nt small RNAs. These results were consistent with previous studies in peanut small RNA sequencing [35, 36, 38].

Fig. 2
figure2

Sequence length distribution of small RNA in different libraries of peanut

More than 85.5% unique small RNAs and 86% total small RNAs were mapped to peanut genome (A. duranensis and A. ipaensis, https://peanutbase.org/) using SOAP2 software (Additional file: Table S2). By blast search against the miRNA database (miRbase 21.0, http://www.mirbase.org/), we identified 113 known miRNAs, belongs to 30 miRNA families (Additional file: Table S3). In addition, clean small RNAs were also used to blast against with GenBank and Rfam database to remove rRNAs, tRNAs, snoRNAs and snRNAs (Additional file: Figure S1). The other small RNAs were used to explore novel miRNAs according the criteria described by the previous studies [35, 39, 40]. In total, we identified 67 novel miRNAs, ranging from 20 nt to 23 nt, and the negative folding free energies of their precursor’s hairpin structures ranged from 184.8 to 19.1 kcal/mol (Additional file: Table S4).

Expression of miRNAs in different peanut varieties response to A. flavus infection

To identify miRNAs in peanut that respond to A. flavus, the differential expressed miRNAs were identified at a significant level (adjusted p-value < 0.05 and fold change ≥1.5) by comparing the normalized miRNA expression level using TP10M (Tags per ten million). Among the 113 known miRNAs identified in peanut, eight and two were up-regulated, twelve and six were down-regulated in response to A. flavus in GT-C20 and Tifrunner, respectively (Fig. 3a, Additional file: Table S3). The total number of differentially expressed known miRNAs in GT-C20 (20) was more than that in Tifrunner (8). Among these differentially expressed known miRNAs, the expression of miR319 were both increased in GT-C20 and Tifrunner after A. flavus inoculation. While the expression of miR157a and miR157d were both decreased in GT-C20 and Tifrunner after treated with A. flavus. In addition, miR157a showed the highest expression abundance in all known miRNAs (Additional file: Table S3). There were also some miRNAs showing different expression trend in response to A. flavus in different peanut varieties. For example, the expression of miR156a was increased in GT-C20 in response to A. flavus, but unchanged in Tifrunner (Additional file: Table S3). Moreover, twelve differentially expressed novel miRNAs were identified in response to A. flavus, including six novel miRNAs from GT-C20 and six novel miRNAs from Tifrunner (Additional file: Table S4). In GT-C20, miRn21, miRn56 and miRn67 were induced or specifically expressed in response to A. flavus, and miRn3, miRn9, miRn18 were depressed in response to A. flavus. In Tifrunner, miRn12, miRn21 and miRn25 were induced or specifically expressed in response to A. flavus, and miRn1, miRn53 and miRn55 were depressed in response to A. flavus (Additional file: Table S4). To confirm the sequencing results, six miRNAs (miR156i, miR166g-3p, miR167a, miR167h, miR172a and miR396e-3p) were selected for stem-loop qRT-PCR analysis (Fig. 4). The expression profiles of these miRNAs were in accordance with the deep sequencing data (Additional file: Figure S2).

Fig. 3
figure3

Numbers of differentially expressed known miRNAs (a) and genes (b) in response to A. flavus

Fig. 4
figure4

Validation of the relative expression level of partial miRNA and mRNA by qRT-PCR. Error bars indicate ± SE obtained from three biological repeats. Student’s T-test was performed to analyze the changes in the gene expression after treated with A. flavus. **denotes the p value < 0.01 and *denotes the p value < 0.05

The miRNA abundance between two different varieties was compared. Before inoculation, the abundance of all 113 known miRNAs of GT-C20 was 245,769 (TP10M), which is lower than that in Tifrunner (615,962). After inoculation, the abundance of all 113 known miRNAs in GT-C20 (269,075.54 TP10M) was also lower than that in Tifrunner (442,962.68 TP10M) (Additional file: Table S3). Under the criterion of adjusted p value (padj) < 0.05 and fold change ≥1.5, a total of 23 and 39 differential expressed known miRNAs between two varieties were identified in control and A. flavus infected samples, respectively (Fig. 3a). Many known miRNAs showed higher expression level in Tifrunner than in GT-C20 in both control and A. flavus infected samples. For examples, in control samples, the expression level of miR157a was 52,711 TP10M in GT-C20 and 291,784 TP10M in Tifrunner. And in treatment samples, the expression level of miR157a-5p was 32,146 TP10M in GT-C20 and 154,326 TP10M in Tifrunner (Additional file: Table S3). In addition, the expression level of miR2118 and miR482a in GT-C20 was also lower than that in Tifrunner (Additional file: Table S3).

Targets identification of miRNAs by high throughout degradome analysis

To validate the targets of peanut miRNAs, a degradome pool was mixed with equal samples from TT, TC, GC and GT. Through high throughput sequencing, 14,338,349 raw reads representing 5,477,891 unique reads were obtained. All reads were used to align against the peanut genome and transcriptome using SOAP2 Program. As a result, 9,580,557 (66.82%) and 6,271,065 (43.73%) reads were successfully mapped in peanut reference genome and transcriptome, respectively (Additional file: Table S5). According to the signature number and abundance of putative cleaved position at each occupied transcript, these cleaved transcripts could be categorized into five classes according to the signature abundance at each occupied transcript position (0, 1, 2, 3 and 4) [41,42,43]. In our dataset, 249, 165, 554, 175 and 402 candidate targets were classified into categories 0, 1, 2, 3 and 4, respectively (Additional file: Table S6). In total, 93 and 176 candidate targets were identified for known and novel miRNAs, respectively (Additional file: Table S7-S8). We found many defense related genes were likely cleaved by miRNAs. For example, two PPR repeat-containing protein genes, Aradu.7H0DM.1 and Aradu.4L72B.1 were potentially cleaved by miR3514 and miRn10, respectively (Additional file: Table S8).

Global mRNA expression profiles in peanut in response to A. flavus growth

To study the target gene expression changes in response to A. flavus growth, eight transcriptome libraries were constructed using the same samples as used in miRNA study. Through transcriptome sequencing, a total of 97,021,101 raw reads were generated with an average of 12 million reads per libraries (Additional file: Table S9). Approximately 77.83% reads were successful mapped to the peanut reference genome, and about 8.18% reads were mapped to multiple regions (Additional file: Table S9). In peanut genome database, 40,636 genes from A. duranensis (A genome) and 46,984 genes from A. ipaensis (B genome) were deposited. In our RNA-seq data, 30,254 (74.45%) genes of A. duranensis and 34,172 (72.73%) genes of A. ipaensis were mapped, indicating that these genes were expressed in peanut seeds.

Differentially expressed genes in response to A. flavus

We further investigated the global gene expression profile of peanut in response to A. flavus. By applying a cut-off criterion of probability ≥0.8 and |log2(fold change) | ≥ 1, a total of 261 (125 up-regulated and 136 down-regulated) and 186 (69 up-regulated and 117 down-regulated) differentially expressed genes (DEGs) were identified in response to A. flavus from Tifrunner and GT-C20, respectively (Fig. 3b, Additional file Figure S3). Among these DEGs in response to A. flavus, 15 DEGs were up-regulated in both GT-C20 and Tifrunner, and 57 DEGs were down-regulated in both GT-C20 and Tifrunner (Fig. 5a). In addition, functional annotation results showed that 139 out of 447 (31.09%) DEGs were defense related genes, such as Chitinase, Heat shock protein, WRKY family transcription factor and temperature-induced lipocalin et al. (Table 1). Moreover, we found that an Indole-3-acetic acid inducible 2 gene (Araip.AN5V8.1) showed opposite expression trend in response to A. flavus. It was up-regulated in resistant variety GT-C20, but down-regulated in susceptible variety Tifrunner (Fig. 5a). The qRT-PCR analysis for several DEGs (Aradu.78 K64.1, Aradu.Q1D8Z.1, Aradu.412P9.1, Aradu.VK4DU.1, Aradu.39 V24.1, and Aradu.DAV01.1) confirmed the RNA-seq results (Fig. 4, Additional file Figure S2).

Fig. 5
figure5

Venn diagram of differentially expressed genes in response to A. flavus (a) and between two peanut varieities (b). Venny analysis was performed using the online software VENNY2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html)

Table 1 Differential expressed defense related genes in response to A. flavus

Differentially expressed genes between resistant and susceptible peanut varieties

We further analyzed the gene expression variation between two peanut varieties. By comparison with control samples (GC and TC), 451 DEGs (346 up-regulated and 105 down-regulated) were identified (Fig. 3b). For treatment samples (GT and TT), 703 DEGs (537 up-regulated and 166 down-regulated) were identified (Fig. 3b). Some disease related genes showed different expression level in two varieties. For example, one NBS-LRR disease resistance gene (Araip.WF303.1, log2 TC/GC = 9.47, log2 TT/GT = 4.66) was highly expressed in Tifrunner, but had a lower expression level in GT-C20. Meanwhile, another NBS-LRR encoding gene (Aradu.168 L7.1) showed the opposite expression trend, highly expressed in GT-C20, but had a lower expression level in Tifrunner (Fig. 3b).

GO and KEGG pathways of DEGs

GO analysis of the DEGs was performed and DEGs could be classified into three categories including molecular function, cellular component and biological process. In category of cellular component, “cell” and “cell part” represented the top terms. For the category of molecular function, the GO term of “catalytic activity” was enriched. For the category of biological processes, the most abundant terms were “metabolic process” and “cellular process” (Fig. 6a). The DEGs between control samples and treatments were also analyzed. Two peanut varieties have difference in many GO terms in both control samples and treatment samples, such as “metabolic process”, “catalytic activity” and “cellular process” (Fig. 6b).

Fig. 6
figure6

GO analysis of differentially expressed genes. a DEGs in response to A. flavus, b DEGs between Tifrunner and GT-C20. Gene are classified into three main categories: biological process, cellular component, and molecular function. The x-axis indicates the number of genes in a category, and the y-axis means the GO tems of genes

To better understand the function and gene regulatory network, KEGG analysis was carried out. The DEGs in response to A. flavus were assigned into 53 and 41 KEGG pathways in Tifrunner and GT-C20, respectively (Fig. 7, Table S1011). The top 16 enriched pathways possibly regulated by A. flavus were summarized (Fig. 7a, Table S10). Most enriched pathways were related to biosynthesis and metabolism, such as “biosynthesis of secondary metabolites”, “flavonoid biosynthesis”, “stilbenoid, diarylheptanoid”, “gingerol biosynthesis”, “metabolic pathways”, and “limonene and pinene degradation” etc. Interestingly, we observed that, in most of the enriched pathways including “biosynthesis of secondary metabolites” and “metabolic pathways”, the numbers of DEGs in Tifrunner was more than that in GT-C20, indicating the difference in metabolic active between two varieties. Interestingly, several pathways of biosynthesis and metabolism were also enriched, such as “biosynthesis of secondary metabolites”, “phenylpropanoid biosynthesis”, “phenylalanine metabolism”, “metabolic pathways”, “flavonoid biosynthesis” and “ascorbate and aldarate metabolism” etc. KEGG results also supported the difference in the metabolic activity between two peanut varieties (Fig. 7b, Table S11).

Fig. 7
figure7

KEGG enrichment analysis of differentially expressed genes. a KEGG enriched pathways in response to A. flavus; b KEGG enriched pathways between Tifrunner and GT-C20

Correlation analysis of miRNAs expression profiles and their target genes

To specify the roles of miRNAs in response to A. flavus growth, the expression profiles of miRNAs and targets were analyzed. In response to A. flavus infection, twenty-two potential miRNAs/targets pairs were identified (Table S12). Degradome sequencing showed these target genes were cleaved by miRNAs. However, the expression profiles of miRNAs were not perfectly negative corresponding to these targets. For example, only five miRNAs/targets pairs showed one-to-one correspondence in Tifrunner, including five up-regulated miRNAs/down-regulated targets and two down-regulated miRNAs/up-regulated targets in response to A. flavus growth. In addition, five miRNAs/targets pairs showed opposite expression trend in Tifrunner after A. flavus infection. For other miRNAs/targets pairs in Tifrunner, only the miRNAs or targets showed differential expression in response to A. flavus, and the one-to-one correspondence expression was not observed.

Between two peanut varieties, forty-four differentially expressed miRNAs/targets pairs were identified (Table S13). We further compared the one-to-one expression correspondence of miRNAs/targets pairs between two peanut varieties. Between control samples and A. flavus infected samples (GC and TC), differentially expressed miRNAs and targets were simultaneously observed only in 16 miRNAs/targets pairs, and these results were not observed in other miRNAs/targets pairs. Among these 16 miRNAs/targets pairs, the expression of miRNA and targets were negatively correlated. For the rest five miRNAs/targets, the expression trend of miRNA and target was in the same trend, both up-regulated or down-regulated. For the flavus infected samples (GT and TT), differentially expressed miRNAs and targets were simultaneously observed in 14 miRNAs/targets pairs, ten of them showed negative correlation in expression of miRNAs and targets.

Discussion

Strategies to reveal the disease resistant mechanism of peanut

Cultivated peanut (Arachis hypogaea. L) is an allotetraploid (2n = 4x = 40, AABB) organism, which emerged as a hybridization of two ancient diploid species, probably A. duranensis (A-genome) and A. ipaensis (B-genome), followed by spontaneous whole genome duplication [44,45,46]. Expressed Sequence Tag (EST) was first used to obtain gene sequence information and their response to biotic and abiotic stresses in peanut [47,48,49]. Based on EST sequencing and Microarray technology, dozens of different genes response to A. flavus infection were identified, which provide a useful information for understanding the molecular mechanism of peanut resistance to A. flavus [1, 6]. Recently, the whole genome sequencing of two diploid wild peanut species has been completed (http://peanutbase.org/). This achievement will greatly promote gene discovery and genetic improvement of this major crop. Here, we employed the high throughput mRNA and miRNA sequencing to investigate the gene expression and regulation of peanut upon A. flavus infection. Our data covered the expression information of more than 64,426 genes. It is the first report of genome-wide gene expression profiling study of peanut in response to A. flavus. Moreover, the integrated analysis of mRNAs and miRNAs might provide insight to understanding the miRNA regulatory in peanut in regulation of the resistance to A. flavus.

Differential expression of defense related genes in resistant and susceptible peanut lines

In the present study, a total of 447 differentially expressed genes (DEGs) were identified in GT-C20 and Tifrunner in response to A. flavus (Fig. 4). Functional annotation results showed that 139 of these genes were defense related genes (Table 1). The majority of defense related DEGs (132 out of 139) were differentially expressed in both peanut varieties, while only seven of the DEGs were differentially expressed only in one peanut variety in response to A. flavus.

Chitin is a key component of fungal cell wall and is considered as an essential factor for maintaining the pathogenicity of the fungus. Chitinases are hydrolytic enzymes that have the activity of degrading chitin. Accumulating evidence indicated that the activity of chitinase was important for plant resistance against pathogens infection [50,51,52], and some chitinases were pathogenesis related (PR) proteins [53,54,55]. We observed that the expression of a chitinase gene (Araip.8C3IU.1) had no change in Tifrunner but was induced in GT-C20. In addition, our data showed that two peroxidases (Aradu.BNR06.1 and Araip.595JK.1) and two copper amine oxidases (Aradu.NCJ0H.1 and Araip.26B5V.1 were up-regulated in GT-C20 but no significant expression difference in Tifrunner (Table 1). Interestingly, peroxidase and blue copper protein genes were also disease related genes [56, 57]. In plant, peroxidases were involved in plant disease resistance response through participating biosynthesis of bactericide and lignin. Copper amine oxidases contribute to catalyze the oxidative de-amination of polyamines and producing hydrogen peroxide (H2O2). In plants, H2O2 derived from polyamine oxidation mediated cell death and was correlated with stress response to pathogen invasion [58, 59]. Plant cell walls were considered as the first security system barrier of defense against bacterial and fungal pathogen. Pectin is an important component in plant cell wall. In plants, pectinesterase/pectinesterase inhibitor acts to modify the cell walls via degration of pectin. In our dataset, two pectinesterase/pectinesterase inhibitor genes (Aradu.S6DQM.1 and Araip.EQZ9W.1) were significantly depressed in Tifrunner but no significant expression difference in response to A. flavus (Table 1). Here, the difference expression trend of these genes in resistant and susceptible peanut varieties provided us a useful clue to reveal the molecular mechanism of peanut response to A. flavus.

Among these 139 defense related DEGs, 18 genes encode cytochrome P450 superfamily genes (Table 1). As a multi-functional gene, the function of cytochrome P450 in defense reaction were revealed in many plants [60, 61]. In rice, over-expressing of CYP71Z2 enhanced resistance to bacterial blight. Further studies confirmed that cytochrome P450 gene contributes to enhance rice disease resistance to pathogens both through mediating diterpenoid phytoalexin accumulation and via suppression of IAA signaling in rice [62].

NBS-LRR genes represent one of the largest disease resistance gene families in plants [63]. So far, the majority of the disease-resistant genes identified by map-based cloning were NBS-LRR genes, such as resistance to rice blast [64], powdery mildew resistant genes Pm3/Pm8 in wheat [65], leaf rust resistance gene Lr10 in wheat [66]. We found two NBS-NBS-LRR type disease resistance genes (Aradu.168 L7.1, Araip.VGW7F.1) were induced in both peanut varieties. However, the increased expression level of the two genes were more dramatic in resistant genotype than that of susceptible genotype peanut. For example, the expression level of Aradu.168 L7.1 was almost equal in both varieties in control samples (log2TC/GC = 0.11). In response to A. flavus, the expression of this gene was induced in both resistant (log2GT/GC = 3.28) and susceptible peanut (log2TT/TC = 0.68) genotype peanuts. After treatment with A. flavus, the expression of this Aradu.168 L7.1 in resistant genotype was significantly higher than that in susceptible genotype (log2GT/TT = 2.72). Araip.VGW7F.1 had the same expression patterns with that of Aradu.168 L7.1, which was induced in both peanut varieties but were more dramatically induced in resistant genotype (log2GT/TT = 2.84). Aradu.168 L7.1 and Araip.VGW7F.1 located in chromosome A05 (1880820–1,893,586) and B05 (1716985–1,720,472) of peanut, respectively, and were considered as orthologous genes. Interestingly, sequence alignment results showed that Aradu.168 L7.1 (Expect: 5e-146, Identities: 316/943(34%) and Positives: 506/943(53%)) and Araip.VGW7F.1 (Expect: 2e-119, Identities: 292/931(31%) and Positives: 490/931(52%)) were all homologous with the disease resistance locus, RPM1, of Arabidopsis. The RPM1 locus conferred resistance to the bacterial pathogen Pseudomonas syringae in Arabidopsis [67, 68] and was supposed to have functional conservation in bean, pea and other crop species [69].

Besides, other several defense related genes also showed differential expression profile between two different peanut genotypes including two cellulose synthase like E1 genes, two WRKY family transcription factor genes, two lipase/lipooxygenase (LOX) genes and one indole-3-acetic acid inducible 2 gene (Table 1). WRKY represents one of largest transcription factor gene families. Some WRKY gene members have been reported in regulating response to various biotic and abiotic stresses. Here, we found that two WRKY genes, Aradu.S7YD6 and Araip.RC4R7, located on chromosome A06 and B06, respectively, were homologous with WRKY33 genes of other plants. In Arabidopsis, WRKY33 was required for resistance to necrotrophic fungal pathogens. In our dataset, the above two WRKY genes were up-regulated in Tifrunner but remained unchanged in GT-C20 which seemed to be the opposite of what was expected in other plants. In addition, we found that two cellulose synthase like E1 genes (Aradu.BZH82, Araip.FL29A) were up-regulated in Tifrunner but did not showed significant differential expression in GT-C20 in response to A. flavus. Cell wall is an important barrier for plants against fungi which also seemed to be the opposite of what was expected in other plants. Previous studies have showed that lipoxygenase pathway was associated with the seed resistance against A. flavus in soybean and maize. For example, inactivation of the ZmLOX3 could increases susceptibility of maize to Aspergillus. Here, two other lipase/lipooxygenase (LOX) genes, Aradu.MAS03 and Araip.081EX were down regulated in the Tifrunner but unchanged in GT-C20. In addition, indole-3-acetic acid inducible genes were involved in response to auxin stimulus and lateral root morphogenesis in other plants. Here, we found that an indole-3-acetic acid inducible 2 gene (Araip.AN5V8.1) was upregulated in the resistant variety and downregulated in the susceptible.

Regulation of miRNAs in peanut response to A. flavus

Due to the sequence similarity, miR2118 was considered as member of miR482 superfamily [34, 70]. Emerging evidence indicated that miR482 and miR2118 both target the P-loop sequence motif of NBS-LRR defense genes in plants, and played a negative regulatory role in plant response to diseases [71]. In soybean, constitutive expression of miR482 led to enhance degradation of targeted R genes, resulting in increased nodulation [72]. In potato, the expression of miR482e was down-regulated in response to Verticillium dahliae infection which resulted in the up-regulation of its NBS-LRR targets. Overexpression of miR482e induced the decrease expression of its NBS-LRR targets, and the transgenic plantlets were highly susceptible to V. dahlia infection [73]. In our datasets, one member of miR2118 and three members of miR482 were identified. We found the overall expression level of miR2118, miR482a, miR482b and miR482b were lower in resistant genotype than in susceptible genotype (Table S3). For example, in resistant genotype, the expression abundance of miR2118 was 74.10 and 91.53 TP10M in mock and treatment samples, respectively. While in susceptible genotype, the expression abundance of miR2118 was 170.22 and 261.62 TP10M in mock and treatment samples, respectively (Table S3). Interestingly, peanut miR482/2118 miRNA family also target NBS-LRR defense genes including Aradu.168 L7, which has the higher expression level in resistant genotype than that in susceptible genotype.

In this study, 66 potential miRNA/targets pairs were identified through a combined analysis of the datasets of small RNA, degradome sequences and mRNA expression profiling. Among them, many SPLs (SQUAMOSA promoter-binding protein-like) might be regulated by miR156/157 miRNA family (Table S1213). For example, Aradu.0GH1S.1 encodes a SPL transcription factor, which is homologous with SPL12 of soybean (LOC100781289) and SPL2 of Arabidopsis (AT5G43270). Degradome sequencing showed that Aradu.0GH1S.1 was potentially cleaved by many miRNA members of miR156/157 family, such as miR156a, miR156e and miR156i, etc. (Table S7). Previous evidences showed that SPL negatively regulated the accumulation of anthocyanin. And further studies proved that miR156/SPL regulated the secondary metabolism and the content of anthocyanin and flavonol in plants [74]. Both anthocyanin and flavonol were important flavonoids. Flavonoids have been shown to have antimicrobial properties against invading microorganisms including A. flavus [75]. In this study, the expression level of Aradu.0GH1S.1 in susceptible genotype was significantly higher than that in resistant genotype, which might be caused by the regulation of miR156. In addition, GO and KEGG results showed that the GO term and pathway related with secondary metabolism of flavonoids was enriched between two peanut varieties. These results provided valuable information for understanding the roles of miR156/SPL in peanut.

Conclusion

Our study reported the integration analysis of the transcriptome, miRNAome and degradome of resistant and susceptible peanut varieties in response to A. flavus. A total of 30 differentially expressed miRNAs, 447 differentially expressed genes and 21 miRNA/targets pairs were identified in response to A. flavus. In addition, a total 62 differentially expressed miRNAs, 451 differentially expressed genes and 44 miRNA/targets pairs were identified between resistant and susceptible peanut varieties. Furthermore, the function of two miRNA/targets regulation pairs were discussed including miR482/2118/NBS-LRR and miR156/157/SPL (Fig. 8). Our study generated a comprehensive dataset for further understand the roles of miRNA in response to A. flavus in peanut.

Fig. 8
figure8

Flow chart for analysis of the transcriptome, miRNAome and degradome of R and S genotype peanut varieties in response to A. flavus. The gray background (top of figure 8) showed the comparison between two peanut varieties. The faint yellow background (bottom of figure 8) showed the comparison between control and treatment by A. flavus

Methods

Plant material and treatment condition

Two cultivated peanuts, including GT-C20 and Tifrunner, were obtained from Shandong Center of Crop Germplasm Resources, Jinan, China, and the voucher specimens were deposited in the Shandong Center of Crop Germplasm Resources. The resistance of GT-C20 and Tifrunner to A. flavus has ever been reported in previous studies [76]. The A. flavus strain (AF-Wh2014) isolated from peanut was kindly provided by lab of Prof. Boshou Liao, Oil Crops Research Institute of Chinese Academy of Agricultural Sciences. The method for evaluating the resistance of the two peanut varieties in response to A. flavus referred to the previous studies and was improved in this study [9, 10, 37]. Briefly, A. flavus was cultured in fresh potato dextrose agar (PDA) medium at 30 °C for 7 d. The conidia were then collected and suspended into spores suspension (106 CFU/ml) with sterile distilled water. For inoculation with A. flavus, peanut seeds were surface sterilized with 70% alcohol for 1 min, and 3% sodium hypochlorite (NaClO) for 10 min. Then, the seed coat was wounded with an area of 6 mm × 5 mm. The wounded area was covered with the same size of a piece of cotton, and 20 μl of spores suspension was added on the cotton. For the control sample, 20 μl sterile water was added. To keep the humidity, the infected seeds were placed in the petri dish with two layers of filter paper wet with 5 ml sterilized water. Finally, seeds were placed into incubator at 30 °C in darkness. Evaluation of peanut resistance to A. flavus was performed as described in a previous study [37]. Infection index was calculated based coverage of A. flavus mycelium in the testa of peanut seed. Peanut interaction with A. flavus was evaluated using the following criteria: high resistant (infection index < 5.0), moderate resistance (infection index ≥5.0 and < 10.0), moderate susceptible (infection index ≥10.0 and < 30.0), susceptible (infection index ≥30.0 and < 50.0).

Libraries construction and deep sequencing

Seeds of GT-C20 and Tifrunner were inoculated with A. flavus for 2 d for small RNA, transcriptome and degradome libraries construction. Un-inoculated seeds were used as control. For each sample, three individually biological replicates were prepared. Total RNA was isolated using Trizol Reagent (Invitrogen, USA) according to the manufacturer’s protocol. RNA quality and integrity was evaluated using electrophoresis on a 1% agarose gel and BioSpectrometer fluorescence (Eppendorf, GER). Then, the total RNAs were used for constructing small RNA, transcriptome and degradome libraries using the methods as described in previous studies, respectively [35, 36, 40, 77]. In brief, for constructing the transcriptome library, the total RNA was purified to obtain the mRNA using oligo (dT) magnetic beads. Then, the mRNA was fragmented into short fragments (about 200 bp) and further used for synthesizing the first strand of cDNA using random hexamer-primer. Next, the second-strand cDNA was synthesized using the buffer, dNTPs, RNase H and DNA polymerase I and then ligated to sequencing adapters. Following the cDNA fragments are enriched by PCR amplification. Finally, the library products are ready for sequencing. For constructing small RNA library, the total RNA was first used for purifying of small RNA molecules (18–30 nt) and then ligated of a pair of adaptors to 5′ and 3′ ends. Following, the sequencing library was prepared through reverse transcription and PCR amplification. For constructing degradome library, following steps were performed. First, mRNA was used as input RNA and annealing with Biotinylated Random Primers. Second, Strapavidin capture of RNA fragments through Biotinylated Random Primers. Third, 5′ adaptor ligation to only those RNAs containing 5-monophosphates. Finally, Reverse transcription and PCR amplification were performed. Before sequencing, the QC step was performed to qualify and quantify of the sample library using Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System. High-throughput sequencing of small RNA, transcriptome and degradome was performed in BGI (Shenzen, China) using Illumina HiSeq2000 platform.

MiRNA identification and bioinformatics analysis

Firstly, raw reads were cleaned up by removing low quality reads, 5′ primer contaminants, reads without 3′ primer, reads with poly A, and reads shorter than 18 nt. The clean reads were aligned by BLAST against Rfam database (version 11.0, http://rfam.janelia.org/) [78, 79] and GenBank noncoding RNA database to remove rRNA, scRNA, snoRNA, snRNA, and tRNA. Then, the rest of reads were aligned with mature miRNAs in miRbase (Version 21, http://www.mirbase.org/) [80, 81] to identify the known miRNAs. Novel miRNAs were identified according to the reported method [39]. The clean reads were first mapped to whole genome sequences of A. duranensis and A. ipaensis (https://peanutbase.org/). Then, the hairpin structures of miRNA precursors were predicted using software Minreap (http://sourceforge.net/projects/mireap) according to the following parameters. The candidate sequences could form a proper secondary hairpin, where mature miRNAs mapped region was with a size of 18–25 nt, and the distance between miRNA and miRNA* is 16–300 nt. The maximal free energy allowed to form a hairpin structure is − 18 kcal/mol [40]. The gene expression level of miRNAs from different sRNA libraries was normalized using TP10M (tags per ten million reads). Differentially expressed miRNAs between different samples were identified using R package DESeq (http://bioconductor.org/packages/DEGseq/) under the criteria of adjusted p value (padj) < 0.05 and fold change ≥1.5 [82].

Degradome sequencing and data analysis

To identify the potential targets regulated by miRNAs, equal amounts of RNAs from control (TC1, TC2, GC1 and GC2) and treatment samples (TT1, TT2, GT1 and GT2) of resistant and susceptible genotypes were mixed together for degradome library construction and deep sequencing. Through preprocessing, clean tags are generated. Then, clean tags were classified by alignment with Genbank, Rfam database, and miRNA database. Next, the reads were mapped to the transcriptome data of A. duranensis and A. ipaensis using SOAP2 program with allowing only two mismatches [83]. The transcriptome data of A. duranensis and A. ipaensis were downloaded from peanutbase (https://peanutbase.org). The sense strand of peanut cDNA was used to predict miRNA cleavage sites using CLeaveLand pipeline [84]. Based on the signature number and abundance of cleaved position at each occupied transcript, the cleaved transcripts could be categorized into five categories (0, 1, 2, 3 and 4) [43].

Transcriptome sequencing and data analysis

To analyze the expression profiles of the target genes, eight independent transcriptome libraries were sequenced. After sequencing, the adaptor sequences, low-quality reads, and empty reads were first removed. Then, the clean reads were aligned with the genome sequences of A. duranensis and A. ipaensis (https://peanutbase.org/) using SOAP2 program [83]. The gene expression level was calculated using FPKM (expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced) method. The relative gene expression level between different samples was calculated using log2 ratio method. Differential expression genes (DEGs) between two samples were identified with the criteria of probability ≥0.8 and |log2(fold change)| ≥ 1 using NOIseq method [85]. To identify the putative biological functions and pathways of the target genes and DEGs, Gene Ontology (GO) annotation and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis were conducted as described previously [77, 86, 87].

Expression validation of miRNA and mRNA using qRT-PCR

To validate the high-throughput sequencing results of miRNAs and target genes, quantitative RT-PCR (qRT-PCR) was performed on ABI 7500 Real-Time PCR System (Applied Biosystems). Three biological replicates were prepared for each sample. For miRNA expression analysis, stem-loop RT-PCR method was used according to the method as described previously [88]. For target genes, primers were designed using primer premier 5.0 software (www.premierbiosoft.com). All primers including the stem-loop miRNA primers, reference genes and targets genes were listed in Additional file: Table S14. The relative expression level between different samples was calculated using the 2-Ct method. Student’s T-test was performed to analyze the changes in gene expression after treated with A. flavus (*P < 0.05, **P < 0.01).

Availability of data and materials

The RNA-seq data in this study were available at NCBI BioProject: PRJNA438019 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA438019).

Abbreviations

GC:

GT-C20 control

GT:

GT-C20 treatment

TC:

Tifrunner control

TT:

Tifrunner treatment

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

NBS-LRR:

Nucleotide-binding site leucine-rich repeat

SPL:

SQUAMOSA promoter binding protein-like

DEG:

Differentially expressed genes

References

  1. 1.

    Guo B, Fedorova ND, Chen X, et al. Gene expression profiling and identification of resistance genes to Aspergillus flavus infection in peanut through EST and microarray strategies. Toxins (Basel). 2011;3:737–53.

  2. 2.

    Guo B, Yu J, Holbrook CC, et al. Strategies in Prevention of Preharvest Aflatoxin Contamination in Peanuts: Aflatoxin Biosynthesis, Genetics and Genomics. Peanut Sci. 2009;36:11–20.

  3. 3.

    Martin CN, Garner RC. Aflatoxin B -oxide generated by chemical or enzymic oxidation of aflatoxin B1 causes guanine substitution in nucleic acids. Nature. 1977;267:863–5.

  4. 4.

    Shen HM, Ong CN, Lee BL, et al. Aflatoxin B1-induced 8-hydroxydeoxyguanosine formation in rat hepatic DNA. Carcinogenesis. 1995;16:419–22.

  5. 5.

    Reddy L, Odhav B, Bhoola K. Aflatoxin B1-induced toxicity in HepG2 cells inhibited by carotenoids: morphology, apoptosis and DNA damage. Biol Chem. 2006;387:87–93.

  6. 6.

    Guo B, Chen X, Dang P, et al. Peanut gene expression profiling in developing seeds at different reproduction stages during Aspergillus parasiticus infection. BMC Dev Biol. 2008;8:12.

  7. 7.

    Wang T, Zhang E, Chen X, et al. Identification of seed proteins associated with resistance to pre-harvested aflatoxin contamination in peanut (Arachis hypogaea L). BMC Plant Biol. 2010;10:267.

  8. 8.

    Wang Z, Yan S, Liu C, et al. Proteomic analysis reveals an aflatoxin-triggered immune response in cotyledons of Arachis hypogaea infected with Aspergillus flavus. J Proteome Res. 2012;11:2739–53.

  9. 9.

    Wang H, Lei Y, Wan L, et al. Comparative transcript profiling of resistant and susceptible peanut post-harvest seeds in response to aflatoxin production by Aspergillus flavus. BMC Plant Biol. 2016;16:54.

  10. 10.

    Wang H, Lei Y, Yan L, et al. Functional Genomic Analysis of Aspergillus flavus Interacting with Resistant and Susceptible Peanut. Toxins (Basel). 2016;8:46.

  11. 11.

    Bertioli DJ, Cannon SB, Froenicke L, et al. The genome sequences of Arachis duranensis and Arachis ipaensis, the diploid ancestors of cultivated peanut. Nat Genet. 2016;48:438–46.

  12. 12.

    Chen X, Li H, Pandey M, et al. Draft genome of the peanut A-genome progenitor (Arachis duranensis) provides insights into geocarpy, oil biosynthesis, and allergens. Proc Natl Acad Sci U S A. 2016;113:6785–90.

  13. 13.

    Reinhart BJ, Weinstein EG, Rhoades MW, et al. MicroRNAs in plants. Genes Dev. 2002;16:1616–26.

  14. 14.

    Jones-Rhoades MW, Bartel DP, Barte B. MicroRNAS and their regulatory roles in plants. Annu Rev Plant Biol. 2006;57:19–53.

  15. 15.

    Chuck G, Candela H, Hake S. Big impacts by small RNAs in plant development. Curr Opin Plant Biol. 2009;12:81–6.

  16. 16.

    Sunkar R, Kapoor A, Zhu JK. Posttranscriptional induction of two cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell. 2006;18:2051–65.

  17. 17.

    Jagadeeswaran G, Saini A, Sunkar R. Biotic and abiotic stress down-regulate miR398 expression in Arabidopsis. Planta. 2009;229:1009–14.

  18. 18.

    Zhang X, Zou Z, Gong P, et al. Over-expression of microRNA169 confers enhanced drought tolerance to tomato. Biotechnol Lett. 2011;33:403–9.

  19. 19.

    Ni Z, Hu Z, Jiang Q, et al. GmNFYA3, a target gene of miR169, is a positive regulator of plant tolerance to drought stress. Plant Mol Biol. 2013;82:113–29.

  20. 20.

    Wang JW, Czech B, Weigel D. miR156-regulated SPL transcription factors define an endogenous flowering pathway in Arabidopsis thaliana. Cell. 2009;138:738–49.

  21. 21.

    Wang JW, Schwab R, Czech B, et al. Dual effects of miR156-targeted SPL genes and CYP78A5/KLUH on plastochron length and organ size in Arabidopsis thaliana. Plant Cell. 2008;20:1231–43.

  22. 22.

    Reyes JL, Chua NH. ABA induction of miR159 controls transcript levels of two MYB factors during Arabidopsis seed germination. Plant J. 2007;49:592–606.

  23. 23.

    Mallory AC, Dugas DV, Bartel DP, et al. MicroRNA regulation of NAC-domain targets is required for proper formation and separation of adjacent embryonic, vegetative, and floral organs. Curr Biol. 2004;14:1035–46.

  24. 24.

    Fang Y, Xie K, Xiong L. Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot. 2014;65:2119–35.

  25. 25.

    Feng H, Duan X, Zhang Q, et al. The target gene of tae-miR164, a novel NAC transcription factor from the NAM subfamily, negatively regulates resistance of wheat to stripe rust. Mol Plant Pathol. 2014;15:284–96.

  26. 26.

    Marin E, Jouannet, Herz A, et al. miR390, Arabidopsis TAS3 tasiRNAs, and their AUXIN RESPONSE FACTOR targets define an autoregulatory network quantitatively regulating lateral root growth. Plant Cell. 2010;22:1104–17.

  27. 27.

    Bowman JL, Class III. HD-zip gene regulation, the golden fleece of ARGONAUTE activity. Bioessays. 2004;26:938–42.

  28. 28.

    Lee DY, An G. Two AP2 family genes, supernumerary bract (SNB) and Osindeterminate spikelet 1 (OsIDS1), synergistically control inflorescence architecture and floral meristem establishment in rice. Plant J. 2012;69:445–61.

  29. 29.

    Wang Y, Wang L, Zou Y, et al. Soybean miR172c targets the repressive AP2 transcription factor NNC1 to activate ENOD40 expression and regulate nodule initiation. Plant Cell. 2014;26:4782–801.

  30. 30.

    Padmanabhan C, Zhang X, Jin H. Host small RNAs are big contributors to plant innate immunity. Curr Opin Plant Biol. 2009;12:465–72.

  31. 31.

    Li F, Pignatta D, Bendix C, et al. MicroRNA regulation of plant innate immune receptors. Proc Natl Acad Sci U S A. 2012;109:1790–5.

  32. 32.

    Boccara M, Sarazin A, Thiebeauld O, et al. The Arabidopsis miR472-RDR6 silencing pathway modulates PAMP- and effector-triggered immunity through the post-transcriptional control of disease resistance genes. PLoS Pathog. 2014;10:e1003883.

  33. 33.

    Li Y, Lu YG, Shi Y, et al. Multiple rice microRNAs are involved in immunity against the blast fungus Magnaporthe oryzae. Plant Physiol. 2014;164:1077–92.

  34. 34.

    Shivaprasad PV, Chen HM, Patel K, et al. A microRNA superfamily regulates nucleotide binding site-leucine-rich repeats and other mRNAs. Plant Cell. 2012;24:859–74.

  35. 35.

    Zhao C, Xia H, Frazier TP, et al. Deep sequencing identifies novel and conserved microRNAs in peanuts (Arachis hypogaea L.). BMC Plant Biol. 2010;10:3.

  36. 36.

    Zhao C, Xia H, Cao T, et al. Small RNA and Degradome deep sequencing reveals Peanut MicroRNA roles in response to pathogen infection. Plant Mol Biol Report. 2015;33:1013–29.

  37. 37.

    Mehan VK, Mcdonald D, Nigam SN, et al. Groundnut cultivars with seed resistant to invasion by Aspergillus flavus. Oleagineux. 1981;36:501–7.

  38. 38.

    Gao C, Wang P, Zhao S, et al. Small RNA profiling and degradome analysis reveal regulation of microRNA in peanut embryogenesis and early pod development. BMC Genomics. 2017;18:220.

  39. 39.

    Meyers BC, Axtell MJ, Bartel B, et al. Criteria for annotation of plant MicroRNAs. Plant Cell. 2008;20:3186–90.

  40. 40.

    Zhang Q, Zhao C, Li M, et al. Genome-wide identification of Thellungiella salsuginea microRNAs with putative roles in the salt stress response. BMC Plant Biol. 2013;13:180.

  41. 41.

    Xu X, Yin L, Ying Q, et al. High-throughput sequencing and degradome analysis identify miRNAs and their targets involved in fruit senescence of Fragaria ananassa. PLoS One. 2013;8:e70959.

  42. 42.

    Yang X, Wang L, Yuan D, et al. Small RNA and degradome sequencing reveal complex miRNA regulation during cotton somatic embryogenesis. J Exp Bot. 2013;64:1521–36.

  43. 43.

    Han X, Yin H, Song X, et al. Integration of small RNAs, degradome and transcriptome sequencing in hyperaccumulator Sedum alfredii uncovers a complex regulatory network and provides insights into cadmium phytoremediation. Plant Biotechnol J. 2016;14:1470–83.

  44. 44.

    Moretzsohn MC, Barbosa AV, Alves-Freitas DM, et al. A linkage map for the B-genome of Arachis (Fabaceae) and its synteny to the A-genome. BMC Plant Biol. 2009;9:40.

  45. 45.

    Moretzsohn MC, Gouvea EG, Inglis PW, et al. A study of the relationships of cultivated peanut (Arachis hypogaea) and its most closely related wild species using intron sequences and microsatellite markers. Ann Bot. 2013;111:113–26.

  46. 46.

    Bertioli DJ, Vidigal B, Nielen S, et al. The repetitive component of the a genome of peanut (Arachis hypogaea) and its role in remodelling intergenic sequence space since its evolutionary divergence from the B genome. Ann Bot. 2013;112:545–59.

  47. 47.

    Bi Y, Liu W, Xia H, et al. EST sequencing and gene expression profiling of cultivated peanut (Arachis hypogaea L). Genome. 2010;53:832–9.

  48. 48.

    Song G, Li M, Xiao H, et al. EST sequencing and SSR marker development from cultivated peanut ( Arachis hypogaea L.). Electron J Biotechnol. 2010;13:7–8.

  49. 49.

    Feng S, Wang X, Zhang X, et al. Peanut (Arachis hypogaea) Expressed Sequence Tag Project: Progress and Application. Comp Funct Genomics. 2012. Article ID 373768. https://doi.org/10.1155/2012/373768.

  50. 50.

    Jiang X, Li Z, Kang Z, et al. Chitinase and plant disease resistance. Acta Agric Boreali-occidentalis Sinica. 2002;11:71–5.

  51. 51.

    van den Burg HA, Westerink N, Francoijs KJ, et al. Natural disulfide bond-disrupted mutants of AVR4 of the tomato pathogen Cladosporium fulvum are sensitive to proteolysis, circumvent Cf-4-mediated resistance, but retain their chitin binding ability. J Biol Chem. 2003;278:27340–6.

  52. 52.

    Cletus J, Balasubramanian V, Vashisht D, et al. Transgenic expression of plant chitinases to enhance disease resistance. Biotechnol Lett. 2013;35:1719–32.

  53. 53.

    Legrand M, Kauffmann S, Geoffroy P, et al. Biological function of pathogenesis-related proteins: four tobacco pathogenesis-related proteins are chitinases. Proc Natl Acad Sci U S A. 1987;84:6750–4.

  54. 54.

    Salzer P, Bonanomi A, Beyer K, et al. Differential expression of eight chitinase genes in Medicago truncatula roots during mycorrhiza formation, nodulation, and pathogen infection. Mol Plant-Microbe Interact. 2000;13:763–77.

  55. 55.

    Gonzalez-Teuber M, Pozo MJ, Muck A, et al. Glucanases and chitinases as causal agents in the protection of Acacia extrafloral nectar from infestation by phytopathogens. Plant Physiol. 2010;152:1705–15.

  56. 56.

    Seevers PM, Daly JM, Catedral FF. The role of peroxidase isozymes in resistance to wheat stem rust disease. Plant Physiol. 1971;48:353–60.

  57. 57.

    Ruan XM, Luo F, Li DD, et al. Cotton BCP genes encoding putative blue copper-binding proteins are functionally expressed in fiber development and involved in response to high-salinity and heavy metal stresses. Physiol Plant. 2011;141:71–83.

  58. 58.

    Rea G, Metoui O, Infantino A, et al. Copper amine oxidase expression in defense responses to wounding and Ascochyta rabiei invasion. Plant Physiol. 2002;128:865–75.

  59. 59.

    Cona A, Rea G, Angelini R, et al. Functions of amine oxidases in plant development and defence. Trends Plant Sci. 2006;11:80–8.

  60. 60.

    Schuler MA, Werck-Reichhart D. Functional genomics of P450s. Annu Rev Plant Biol. 2003;54:629–67.

  61. 61.

    Mao G, Seebeck T, Schrenker D, et al. CYP709B3, a cytochrome P450 monooxygenase gene involved in salt tolerance in Arabidopsis thaliana. BMC Plant Biol. 2013;13:169.

  62. 62.

    Li W, Wang F, Wang J, et al. Overexpressing CYP71Z2 enhances resistance to bacterial blight by suppressing auxin biosynthesis in rice. PLoS One. 2015;10:e0119867.

  63. 63.

    Torii KU. Leucine-rich repeat receptor kinases in plants: structure, function, and signal transduction pathways. Int Rev Cytol. 2004;234:1–46.

  64. 64.

    Qu S, Liu G, Zhou B, et al. The broad-spectrum blast resistance gene Pi9 encodes a nucleotide-binding site-leucine-rich repeat protein and is a member of a multigene family in rice. Genetics. 2006;172:1901–14.

  65. 65.

    Hurni S, Brunner S, Buchmann G, et al. Rye Pm8 and wheat Pm3 are orthologous genes and show evolutionary conservation of resistance function against powdery mildew. Plant J. 2013;76:957–69.

  66. 66.

    Feuillet C, Travella S, Stein N, et al. Map-based isolation of the leaf rust disease resistance gene Lr10 from the hexaploid wheat (Triticum aestivum L.) genome. Proc Natl Acad Sci U S A. 2003;100:15253–8.

  67. 67.

    Stahl EA, Dwyer G, Mauricio R, et al. Dynamics of disease resistance polymorphism at the Rpm1 locus of Arabidopsis. Nature. 1999;400:667–71.

  68. 68.

    Mackey D, Holt BR, Wiig A, et al. RIN4 interacts with Pseudomonas syringae type III effector molecules and is required for RPM1-mediated resistance in Arabidopsis. Cell. 2002;108:743–54.

  69. 69.

    Dangl JL, Ritter C, Gibbon M, et al. Functional homologs of the Arabidopsis RPM1 disease resistance gene in bean and pea. Plant Cell. 1992;4:1359–69.

  70. 70.

    Zhu QH, Fan L, Liu Y, et al. miR482 regulation of NBS-LRR defense genes during fungal pathogen infection in cotton. PLoS One. 2013;8:e84390.

  71. 71.

    de Vries S, Kloesges T, Rose LE. Evolutionarily dynamic, but robust, targeting of resistance genes by the miR482/2118 gene family in the Solanaceae. Genome Biol Evol. 2015;7:3307–21.

  72. 72.

    Li H, Deng Y, Wu T, et al. Misexpression of miR482, miR1512, and miR1515 increases soybean nodulation. Plant Physiol. 2010;153:1759–70.

  73. 73.

    Yang L, Mu X, Liu C, et al. Overexpression of potato miR482e enhanced plant sensitivity to Verticillium dahliae infection. J Integr Plant Biol. 2015;57:1078–88.

  74. 74.

    Gou JY, Felippes FF, Liu CJ, et al. Negative regulation of anthocyanin biosynthesis in Arabidopsis by a miR156-targeted SPL transcription factor. Plant Cell. 2011;23:1512–22.

  75. 75.

    Medina ML, Kiernan UA, Francisco WA. Proteomic analysis of rutin-induced secreted proteins from Aspergillus flavus. Fungal Genet Biol. 2004;41:327–35.

  76. 76.

    Zhang H, Scharfenstein L, Carter-Wientjes C, et al. Lack of aflatoxin production by Aspergillus flavus is associated with reduced fungal growth and delayed expression of aflatoxin pathway genes. World Mycotoxin J. 2014;8:335–40.

  77. 77.

    Wang J, Zhang Q, Cui F, et al. Genome-wide analysis of gene expression provides new insights into cold responses in Thellungiella salsuginea. Front Plant Sci. 2017;8:713.

  78. 78.

    Griffiths-Jones S, Moxon S, Marshall M, Khanna A, Eddy SR, Bateman A. Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res. 2005;33(suppl_1):D121–4.

  79. 79.

    Burge SW, Daub J, Eberhardt R, Tate J, Barquist L, Nawrocki EP, Eddy SR, Gardner PP, Bateman A. Rfam 11.0: 10 years of RNA families. Nucleic Acids Res. 2012;41(D1):D226–32.

  80. 80.

    Griffiths-Jones S. The microRNA registry. Nucleic Acids Res. 2004;32(suppl_1):D109–11.

  81. 81.

    Griffiths-Jones S, Grocock RJ, Van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(suppl_1):D140–4.

  82. 82.

    Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2009;26(1):136–8.

  83. 83.

    Li R, Yu C, Li, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25:1966–7.

  84. 84.

    Addo-Quaye C, Miller W, Axtell MJ. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2008;25(1):130–1.

  85. 85.

    Tarazona S, Garcíaalcalde F, Dopazo J, et al. Differential expression in RNA-seq: a matter of depth. Genome Res. 2011;21:2213–23.

  86. 86.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25.

  87. 87.

    Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

  88. 88.

    Varkonyi-Gasic E, Wu R, Wood M, et al. Protocol: a highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods. 2007;3:12.

Download references

Acknowledgments

The authors thank Prof. Boshou Liao and Prof. Yong Lei (Oil Crops Research Institute of Chinese Academy of Agricultural Sciences) for providing the A. flavus strain (AF-Wh2014).

Funding

This work was supported by the Natural Science Foundation of China [31861143009]; Shandong Province Germplasm Innovation and Utilization Project [2016LZGC025], Shandong Province Natural Science Foundation [ZR2017MC005, ZR2015YL061, ZR2015CQ012 and ZR2018BC029], Shandong Provincial Key Research and Development Program (2019GNC106026 and 2019GSF107008) and Taishan Scholar Project of Shandong Province (ts20190964); Agricultural Scientific and Technological Innovation Project [CXGC2018E13, CXGC2016B02] and Young Talents Training Program of Shandong Academy of Agricultural Sciences. The founders did not play any roles in the design, analysis, and interpretation of this study or relevant data.

Author information

Affiliations

Authors

Contributions

CZ and XW conceived and designed the experiment; CZ, TL, YZ1, BZ, JQ, SF and BG analyzed the data and wrote the manuscript; TL, AL, SZ, LH, HX, PL and YZ2 performed the experiment. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Xingjun Wang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

All authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1. Distinct RNA fragment categories in each library.

12870_2020_2426_MOESM2_ESM.pptx

Additional file 2: Figure S2. Comparison of sequencing and qRT-PCR results for the miRNAs and genes using heatmap. Heatmap was generated by online software Morpheus (https://software.broadinstitute.org/morpheus) according to the relative expression in comparison with the control sample TC.

Additional file 3: Figure S3. Volcano plots showing the gene expression differences among different samples. The RPKM-normalized transcript count data sets were analysed by NOIseq: the x-axis shows the log-ratio (gene expression fold change after challenge) the y-axis the probability for each gene of being differentially expressed. The up-regulated and down-regulated genes were identified under the criterion of probability ≥0.8 and |log2(fold change)| ≥ 1.

Additional file 4: Table S1. Summary of small RNA reads from the individual libraries.

Additional file 5: Table S2. Statistical analysis of total sRNA mapped in peanut genome.

Additional file 6: Table S3. Known miRNA information in two peanut varieties.

Additional file 7: Table S4. Information of novel miRNAs identified from peanut miRNA libraries.

Additional file 8: Table S5. Classification analysis of reads from degradome library.

Additional file 9: Table S6. Categories of candidate cleaved sites.

Additional file 10: Table S7. Target genes of known miRNAs identification by degradome sequencing.

Additional file 11: Table S8. Target genes of novel miRNAs identification by degradome sequencing.

Additional file 12: Table S9. Summary of RNA-seq reads from the individual libraries.

Additional file 13: Table S10. Top sixteen KEGG pathways in Tifrunner and GT-C20 in response to AF infection.

Additional file 14: Table S11. Top sixteen enriched KEGG pathways between Tifrunner and GT-C20.

Additional file 15: Table S12. Correlation analysis of miRNAs expression profiles and their target genes in response to A. flavus.

Additional file 16: Table S13. Correlation analysis of miRNAs expression profiles and their target genes between peanut varieties.

Additional file 17: Table S14. Primers used in this study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Zhao, C., Li, T., Zhao, Y. et al. Integrated small RNA and mRNA expression profiles reveal miRNAs and their target genes in response to Aspergillus flavus growth in peanut seeds. BMC Plant Biol 20, 215 (2020). https://doi.org/10.1186/s12870-020-02426-z

Download citation

Keywords

  • Peanut
  • Aspergillus flavus
  • microRNA
  • Transcriptome
  • Degradome