MicroRNA Expression Profile during Aphid Feeding in Chrysanthemum (Chrysanthemum morifolium)

MicroRNAs (miRNAs) are important regulators of gene expression, affecting many biological processes. As yet, their roles in the response of chrysanthemum to aphid feeding have not been explored. Here, the identity and abundance of miRNAs induced by aphid infestation have been obtained using high-throughput Illumina sequencing platform. Three leaf small RNA libraries were generated, one from plants infested with the aphid Macrosiphoniella sanbourni (library A), one from plants with mock puncture treatment (library M), and the third from untreated control plants (library CK). A total of 7,944,797, 7,605,251 and 9,244,002 clean unique reads, ranging from 18 to 30 nucleotides (nt) in length, were obtained from library CK, A and M, respectively. As a result, 303 conserved miRNAs belonging to 276 miRNAs families and 234 potential novel miRNAs were detected in chrysanthemum leaf, out of which 80, 100 and 79 significantly differentially expressed miRNAs were identified in the comparison of CK-VS-A, CK-VS-M and M-VS-A, respectively. Several of the differentially abundant miRNAs (in particular miR159a, miR160a, miR393a) may be associated with the plant's response to aphid infestation.


Introduction
Small RNAs (smRNAs), ranging in length from 18 to 40 nucleotides (nt), regulate gene expression via their interaction with specific mRNA targets, thereby playing a significant role throughout plant growth and development [1,2]. MicroRNAs (miRNAs), approximately 22 nt in size, are the most abundant species of endogenous noncoding smRNAs in plants [3,4]. They have been implicated in the response to both biotic and abiotic stress [5,6], by guiding the cleavage or repressing the translation of their target mRNAs, which is based on nearly perfect complementarity with the targets, resulting in gene silencing at post-transcriptional levels [3,7].
The detection of miRNAs brings us into a new vision in the comprehension of gene regulation in plants. Evidence from the discovery that miR398 targeting two Cu/Zn superoxide dismutases (CSD1 and CSD2) and a copper chaperone for CSD1 responds to both oxidative stress quality reads and those harboring adaptor sequence, a total of 7,944,797, 7,605,251 and 9,244,002 clean unique reads, ranging from 18 to 30 nucleotides (nt) in length, were obtained from library CK, A and M, respectively (Table 1). All clean reads were deposited in the NCBI Sequence Read Archive database (http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?) under accession number SRP051318. (Accession number SRS796738 for library CK; Accession number SRS799803 for library A; Accession number SRS799804 for library M). The proportion of clean reads was >99.49% in each library. As seen in Fig 1, the most abundant length was 24 nt, just as for rice, Arabidopsis and soybean. The three libraries each harbored a similar proportion of miRNAs (about 0.281% of unique reads and 5.021% of the total reads), rRNA (1.437% and 7.765%), snRNA (0.014% and 0.009%), snoRNA (0.005% and 0.004%) and tRNA (0.175% and 1.516%). The remaining unannotated reads were used to identify the new miRNAs.

Identification of conserved miRNAs and their expression profiles
By matched to known miRNAs from miRBase database, a total of 303 conserved miRNAs belonging to 276 miRNAs families were identified in the three libraries (S1 Table). Of these, 222 conserved miRNAs were detected in library CK, and 221 and 224 miRNAs in library A and M, respectively. Among them, 34, 24 and 34 miRNAs were specifically expressed in library CK, A and M, respectively; 174, 176 and 167 miRNAs were co-expressed in library CK and A, library A and M or library CK and M, respectively; and 153 miRNAs were simultaneously expressed in library CK, A and M (Fig 2A, S1 Table).

Prediction of conserved miRNA targets and mapping of cleavage sites in the targets
Based on the characteristic of near-perfect complementarity with their targets of plant miR-NAs, the transcriptome of chrysanthemum was used to predict targets of the conserved miR-NAs according to Allen et al. [24,25]. In the comparison between CK and A (CK-VS-A), 26 out of 80 differentially expressed miRNAs were predicted to target 145 genes (26/80/145). And the corresponding numbers in CK-VS-M and M-VS-A were 38/100/159 and 27/79/110, respectively (S3 Table). For examining the influences of the conserved miRNAs on their predicted target genes, RNA ligase-mediated 5' rapid amplification of cDNA ends (RLM-5' RACE) was performed to map internal cleavage sites in the targets using a GeneRacer kit (Invitrogen, Carlsbad, CA). The miR160a-guided cleavage of Unigene29593_All which encodes an ARF gene and the miR393a-guided cleavage of CL1261.Contig1_All encoding one TRANS-PORT INHIBITOR RESPONSE 1 (TIR1) protein were validated in chrysanthemum (Fig 4).

Gene ontology (GO) functional classification
To investigate the potential roles of miRNAs, we performed GO analysis in terms of the targets of conserved miRNAs that were differentially expressed. In CK-VS-A comparison, 88 target genes were categorized as "biological process", 62 as "cellular component" and 43 as "molecular function" (88/62/43), and 97/70/59 and 90/83/40 in CK-VS-M and M-VS-A, respectively ( Fig 5).
As shown in Fig 5, most of target genes were categorized into "biological process", and majority of them were associated with cellular process, metabolic process and response to stimulus in CK-VS-A, CK-VS-M and M-VS-A. Furthermore, in terms of cellular component, the majorities were involved in cell, cell part and organelle, and most of the targets were involved in binding and catalytic activity in terms of molecular function. Moreover, on the aspect of numbers of target genes in the categories, these genes were mainly enriched in binding and catalytic activity in CK-VS-A and CK-VS-M. Prediction of potential novel miRNAs In total, 234 potential novel miRNAs were identified in the three libraries (S4 Table). Of these, 137 potential novel miRNAs were detected in library CK, and 123 and 135 miRNAs in library A and M, respectively. Among them, 49, 37 and 46 miRNAs were specifically expressed in library CK, A and M, respectively; 72, 73 and 75 miRNAs were co-expressed in library CK and A, library A and M or library CK and M, respectively; and 59 miRNAs were simultaneously expressed in library CK, A and M ( Fig 6A, S4 Table). The precursor sequences of these potential novel miRNAs were listed in S1 File. According to fold-change (|log 2 Ratio(treatment/control)| ! 1.0) and P-value (< 0.05), there were 30 significantly differentially expressed miRNAs (11 miRNAs up-regulated and 19 miR-NAs down-regulated, 11/19) in library CK and A comparison (CK-VS-A), and 26 (8/18) and 20 (11/9) miRNAs in CK-VS-M and M-VS-A, respectively (S5 Table) Table).

miRNAs responding to aphid and mock puncture treatments
Of the 537 miRNAs, 153 conserved miRNAs and 59 potential novel miRNAs were found in all three treatments, occupying 39.5% of the total (Figs 2A and 6A). 24 conserved miRNAs and 37 potential novel miRNAs were specifically detected in plants subjected to aphid herbivory. There were 80 puncture-specific miRNAs, in which 34 miRNAs were conserved and 46 miR-NAs were novel. 23 conserved miRNAs and 14 potential novel miRNAs were specific to both aphid and mock puncture treatments. Compared with the control (CK) (Figs 2B and 6B), 32 conserved miRNAs and 13 potential novel miRNAs were differentially responded to aphid feeding. 48 conserved miRNAs and 17 potential novel miRNAs differentially response to both aphid infestation and mock puncture were identified. Additionally, 52 conserved miRNAs and 9 potential novel miRNAs were differentially expressed during mock puncture treatment.  To validate the results of Illumina deep sequencing, several miRNAs from library CK and A (CK: control; A: aphid infestation treatment) were chosen randomly for stem-loop qRT-PCR. For comparison of fold change between Illumina sequencing and qRT-PCR, scatterplots were generated using the log 2 fold change determined by Illumina sequencing and qRT-PCR. As shown in Fig 7, the qRT-PCR results revealed that the expression of these miRNAs showed a similar tendency (R 2 = 0.57) [26] with the Illumina sequencing data, suggesting the reproducibility and accuracy of sequencing results.

Discussion
As important regulators of gene expression, miRNAs are widely discovered in animals and plants. Thousands of miRNAs of plants were annotated in the miRBase database, including rice and Arabidopsis [27,28]. Here, the capacity of the Illumina sequencing platform to generate large quantities of sequence in a single run has been exploited to achieve a global picture of the miRNA content of chrysanthemum, a larger genome polyploid, and in particular to identify those which may be involved in the plant's response to aphid infestation. As a result, a total of 303 conserved miRNAs belonging to 276 miRNAs families and 234 potential novel miRNAs were identified in three libraries constructed from pooled leaf tissues of Chrysanthemum morifolium 'nannongxunzhang' that were collected at different time points with (A) or without (CK) aphid infestations and mock puncture treatment (M) (S1 and S4 Tables). In this study, the abundance of a large number of miRNAs was increased in response both to aphid feeding (CK-VS-A) and mock puncture treatments (CK-VS-M). To figure out the potential role of these miRNAs, several specific miRNAs that are candidates for the resistance to aphids were further analyzed. Chrysanthemum miR159a induced by aphid feeding and mock puncture treatments was predicted to target a GAMYB-like 2 gene. The GAMYB or GAMYB-like genes encode R2R3 MYB domain transcription factors involving in gibberellin (GA) signaling pathway [29]. In Arabidopsis, miR159 suppresses the expression of GAMYBlike gene(s) regulating cell proliferation in vegetative tissues and promoting programmed cell death in the aleurone layer [30]. Pegadaraju and co-workers [31] reported that the feeding of the aphid Myzus persicae on A. thaliana up-regulated the expression of a senescence-associated gene (SAG), triggering the premature senescence of the aphid-infested leaves, thereby restricting the multiplication of the aphids, contributing to the defense against M. persicae in A. thaliana. Aphid feeding could create a strong sink in insect-infested tissue, leading to increased flow of nutrients to those tissues in the host [32]. So we propose that miR159a regulation might be implicated in the PCD induced by aphid infestation (aphid saliva and mechanical injury), which redirects the resource allocation and contributes to resist aphids.  Both miR160a and miR393a were down-regulated by both aphid infestation and the mock puncture treatment. The RLM 5'-RACE assay successfully identified their targets in chrysanthemum as members of Auxin Response Factors (ARFs) family and TRANSPORT INHIBITOR RESPONSE 1 (TIR1) (Fig 4), respectively. The transcript abundance of the two targets was both increased in response to aphid feeding (S1 and S2 Figs). TIR1 is part of the ubiquitin-ligase complex SCF TIR1 mediating the degradation of Aux/IAA proteins that repress auxin signaling through heterodimerization with ARFs [33,34]. Navarro et al. [10] showed that miR393 induced by a flagellin-derived peptide negatively regulated mRNAs for the F-box auxin receptors TIR1, AFB2 and AFB3, thereby repressing the auxin signaling and restricting Pseudomonas syringae growth, indicating auxin in the susceptibility of disease and miRNA-mediated suppression of auxin signaling in resistance. In different eco-types of Arabidopsis, Kuśnierczyk and co-workers [35] found that several genes which belong to the auxin biosynthesis pathway were up-regulated after green peach aphid (M. persicae) infestation. The induced down-regulation of miR160a and miR393a may suggest the involvement of miRNA-mediated gene regulation in the chrysanthemum-aphid interaction.
Unfortunately, among those conserved miRNAs, only half of them could be predicted to targets, of which two were further validated through RLM-5' RACE method. The reasons that only a small number of validated targets may be as follows: the limited coverage of the chrysanthemum transcriptome; the limitation of the RLM-5' RACE method which is hard to detect the low abundance of the cleavage products and is difficult to identify genes that are targeted by more than one miRNA; and in some cases, miRNAs might silence gene by repressing their translation rather than their transcription. Although there is a deal of conservation across species in terms of miRNA targets, it has been noted that some known miRNAs target additional, non-conserved genes in different species [4,36]. An example is provided by miR159a, which in addition to GAMYB-like 2, also interacts with a gene encoding a 50S ribosomal protein. Similarly, miR397a targets a glucosidase gene in addition to TIR1. Degradome sequencing [37] and transient coexpression of miRNAs and their predicted targets in Nicotiana benthamiana leaves [38] both helps the prediction of target genes, which will be our ongoing study. Moreover, with the development of sequencing technologies, the completion and release of chrysanthemum genome and transcriptome sequence will be expectable and would aid in the identification of predicted conserved and additional target genes.
Furthermore, the mock puncture treatment was included in the present experiments in an attempt to separate the plant's wounding response caused by the insertion of the aphid stylet from its response to feeding. In spite of some differentiations between aphid stylet and puncture, for example, aphid stylets were often wrapped by saliva that contains a complicated mixture of enzymes and can lead to defense responses [39]. Also, the mechanical degree of puncture treatment should be different from aphid stylets. Results show that it does have some similarities between aphid feeding and puncture treatment. As seen in Figs 2B and 6B, 32 conserved miRNAs and 13 potential novel miRNAs were specifically differentially expressed in CK-VS-A; whereas 48 conserved miRNAs and 17 potential novel miRNAs were co-expressed in CK-VS-A and CK-VS-M, suggesting that miRNAs co-expressed in response to aphid feeding and puncture treatment might be involved in wound-induced responses by aphid, otherwise miRNAs may specifically respond to aphid sucking. These will allow us to find out the potential effects of aphid stylets and refine the processes of defense responses.
Here, we mainly focus on the discussions of conserved miRNAs. In the absence of whole genome annotation, the chrysanthemum transcriptome was also used to predict targets for potential novel miRNAs according to the criteria used by Allen et al. [24,25]. However, the sequences of potential targets of miRNAs were too short to get reliable prediction.
Our results indicated that the interaction between chrysanthemum and Macrosiphoniella sanbourni shared some common with interactions between other plants and aphids in terms of changes in expression of miRNAs, such as miR159a, miR160a and miR393a which have been previously proposed to be involved in response to aphid attack in Arabidopsis and melon [4,40], suggesting that present study could provide valuable clues for dissecting other plants' responses to phloem feeding insects. We speculate that these miRNAs might be associated with the aphid infestation in chrysanthemum based on their expression profiles, and the further validations of the two targets of miR160a and miR393a in chrysanthemum suggested that these miRNAs are of worth to be studied further. Additional evidence should be collected to verify if these miRNAs contribute to the resistance of aphid in chrysanthemum, which will be our ongoing research.

Conclusions
In summary, this study has demonstrated that both aphid infestation and mock puncture treatment induced a large change in the spectrum of miRNAs present in the chrysanthemum leaf. Some of the induced miRNAs are likely to be associated with the plant's defense response to aphid feeding, and thus are of interest in efforts to genetically manipulate aphid resistance in chrysanthemum. The number of differentially abundant miRNAs uncovered implies that multiple defense pathways are likely activated by aphid infestation, allowing the plant a degree of genetic flexibility in combating this form of biotic stress.

Plants growth
Cuttings of the aphid resistant Chrysanthemum morifolium cultivar 'Nannongxunzhang' were obtained from the Chrysanthemum Germplasm Resource Preserving Centre, Nanjing Agricultural University, China. Seedlings were grown in a 2:1 mixture of garden soil and vermiculite without added fertilizer under 80% relative humidity, a 16 h photoperiod (160 μmol m -2 s -1 photon flux density), and a day/night temperature of 25/18°C in a greenhouse. Uniform plants at the 6-8 leaf stage were used in all experiments.

Aphid infestation and acupuncture treatment
Aphids (Macrosiphoniella sanbourni Gillette) were collected from field-grown chrysanthemum plants and second instar nymphs were raised for the infestation experiment (A). Because the aphids were from our own chrysanthemum planting base in 'Suoshi' village of nanjing in China, no specific permissions were required for any locations/activities. Aphid infestation (A) and mock puncture treatment (M) were carried according to Xia et al. [41]. For aphid infestation treatment (A), a set of 20 aphids was placed on the third fully expanded leaf (counting from the top of the stem), which was then enclosed by transparent, ventilated plastic cylinder (2 cm height × 5 cm diameter) and sealed at the base of the petiole. The same enclosure was also used for the control (CK) and mock puncture treatment (M). For the mock puncture treatment (M), designed to partially simulate the mechanical stress resulting from aphid penetration, the third fully expanded leaf of each plant was punctured 5 times at 0 h, 10 times at 24 h, and 15 times at 48 h with a needle (approximately 0.30 mm diameter) [42]. Leaves of three individual seedlings for each treatment were harvested at 0 h, 3 h, 6 h, 12 h, 24 h and 48 h. Before harvest, aphids were removed by spraying with 1% (v/v) SDS solution, which caused aphids to remove their mouthparts from plant tissues and then removed the aphids from the leaves by flushing the plants with deionized water. Harvested materials were immediately frozen in liquid nitrogen and stored at -80°C for the following experiments. The samples collected at defined time points of each treatment were pooled for sequencing.
RNA extraction, small RNA library construction and Illumina deep sequencing Total RNA from leaf tissue of three separate libraries (CK, M, A) was extracted using RNAiso reagent (TaKaRa, Otsu, Japan), following the manufacturer's instructions. The integrity and quality of the total RNA was evaluated using a 2100 bioanalyzer RNA Nano chip device (Agilent, Santa Clara, CA, USA) and agarose gel electrophoresis, and the concentration was measured with a ND -1000 spectrophotometer (NanoDrop, Wilmington, DE).
After the total RNA isolation, low molecular weight (LMW) RNAs (18-30 nt) were isolated by polyacrylamide gel electrophoresis. The LMW RNAs were then ligated to a 5' RNA adapter and a 3' RNA adapter. A reverse transcription reaction followed by low cycle PCR amplification (initial denaturation at 98°C for 30 s, followed by 15 cycles of 98°C for 10 s, 60°C for 30 s and 72°C for 15 s, and then 72°C for 10 min) was performed to generate sufficient template for SBS-sequencing by synthesis, which can decrease the loss of nucleotides caused by the secondary structure. PCR products approximately 100 bp were collected by gel purification and sequenced on Illumina Genome Analyzer II platform according to manufacturer's instructions (Illumina Inc, USA).

Bioinformatic analysis of sequencing data
The raw reads from Illumina Genome Analyzer II were initially processed to remove adaptor sequences, low-quality reads as well as contaminants to get clean reads. Reads ranging from 18-25 nt in length were matched to known miRNAs from miRBase database (http://www. mirbase.org/) to identify conserved miRNA. To remove reads from ncRNA (noncoding RNA), such as rRNAs (ribosomal RNAs), tRNAs (transfer RNAs), snRNAs (small nuclear RNAs), snoRNAs (small nucleolar RNAs) and repeat RNA, we aligned the reads with Rfam (http:// www.sanger.ac.uk/) and Genbank database (http://www.ncbi.nlm.nih.gov/). Sequences matching exons and introns of mRNA were also removed to avoid mRNA contamination using overlap software developed by BGI TechSolutions (Shenzhen, China). The remaining reads without any annotation were used to predict potential novel miRNAs by Mireap (http://sourceforge. net/projects/mireap) according to Friedländer et al. [43]. And because of the lack of genome sequence, EST and transcriptome sequences were used to find potential novel miRNAs in chrysanthemum. In the absence of whole genome annotation, the transcriptome of chrysanthemum was used to predict targets of the identified miRNAs according to the criteria used by Allen et al. [24,25]. All sequencing data have been deposited at the sequence read archive (SRA) of NCBI.
The expression abundance of each miRNA was estimated by std, normalized expression level of miRNA in a sample. Normalization formula: Normalized expression = (Original counts of miR-NAs / Total count of clean reads) Ã 10 6 . If the original miRNA expression in a library was zero, the normalized expression was adjusted to 0.01; if the miRNA gene expression of the two samples was less than 1, it was not used to analyze the differential expression due to its low expression. The fold change between two samples was calculated as: fold change = log 2 (treatment/control). The comparisons of the miRNA expression levels between the three libraries were performed using Poisson distribution (P < 0.05). P-value formula: (N 1 : the total clean tag number of sample 1; N 2 : total clean tag number of sample 2; x: tag number of gene A in sample 1; y: tag number of gene A in sample 2). The criteria used for assigning significance were: P-value < 0.05 and estimated absolute |log 2 Ratio(treatment/control)| ! 1.0. Three replicates were used for analysis. To clearly visualize the expression profile data, a heatmap was generated using Cluster v3.0 and Treeview software.

Validation of the expression of miRNAs by stem-loop qRT-PCR
A well-developed stem-loop quantitative real-time PCR (qRT-PCR), which could distinguish single nucleotide deference between miRNAs, was performed to validate the expression of miRNAs identified by sequencing. It includes two steps: a reverse transcription reaction and a real time quantitative PCR reaction. Reverse transcription reactions were carried according to Chen et al [44]. The samples were collected as described above. All primers involved are listed in S6 Table. PCR was performed using a Eppendorf AG 22331 Hamburg thermocycler in a 20 ul volume containing 10 ul SYBR Green PCR master mix (TaKaRa, Japan), 0.2 uM of each primer (S6 Table) and 10 ng cDNA, and the amplification programme including an initial denaturation at 95°C for 60 s, followed by 40 cycles of 95°C for 15 s, 60°C for 15 s and 72°C for 20 s). At the end of the cycling process, a melting-curve analysis from 55 to 95°C with a heating rate of 0.5°C s -1 was performed to determine specificity of amplified products. The chrysanthemum EF1α gene was used as a reference. Relative expression levels were calculated using the 2 -ΔΔCT method [45]. Three independent biological replicates of each sample and three technical replicates of each biological replicate were used for qRT-PCR analysis.

RLM-5' RACE
RNA ligase-mediated 5' rapid amplification of cDNA ends (RLM-5' RACE) was carried to map cleavage sites in the predicted target genes of miRNAs according to Llave et al. [46] using a Gen-eRacer kit (Invitrogen, Carlsbad, CA). mRNA was purified from total RNA extracted from leaf tissues using PolyATtract mRNA Isolation System IV kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. The final PCR products were gel purified and connected to pMD19-T (TaKaRa, Japan) for sequencing. At least 15 independent clones were sequenced.

Accession numbers
The data sets supporting the results of this article are available in the NCBI Sequence Read Archive (SRA) database under accession number SRP051318, http://www.ncbi.nlm.nih.gov/ sra/?term=SRP051318. The expression abundance of each miRNA was estimated by std, normalized expression level of miRNA in a sample. If the original miRNA expression in a library was zero, the normalized expression was adjusted to 0.01. The criteria used for assigning significance were: P-value < 0.05 and estimated absolute |log2Ratio(treatment/control)| ! 1.0. "-": no statistical significance. (XLS) S6