Coordinated miRNA/mRNA Expression Profiles for Understanding Breed-Specific Metabolic Characters of Liver between Erhualian and Large White Pigs

MicroRNAs (miRNAs) are involved in the regulation of various metabolic processes in the liver, yet little is known on the breed-specific expression profiles of miRNAs in coordination with those of mRNAs. Here we used two breeds of male newborn piglets with distinct metabolic characteristics, Large White (LW) and Erhualian (EHL), to delineate the hepatic expression profiles of mRNA with microarray and miRNAs with both deep sequencing and microarray, and to analyze the functional relevance of integrated miRNA and mRNA expression in relation to the physiological and biochemical parameters. EHL had significantly lower body weight and liver weight at birth, but showed elevated serum levels of total cholesterol (TCH), high-density lipoprotein cholesterol (HDLC) and low-density lipoprotein cholesterol (LDLC), as well as higher liver content of cholesterol. Higher serum cortisol and lower serum insulin and leptin were also observed in EHL piglets. Compared to LW, 30 up-regulated and 18 down-regulated miRNAs were identified in the liver of EHL, together with 298 up-regulated and 510 down-regulated mRNAs (FDR<10%). RT-PCR validation of some differentially expressed miRNAs (DEMs) further confirmed the high-throughput data analysis. Using a target prediction algorithm, we found significant correlation between the up-regulated miRNAs and down-regulated mRNAs. Moreover, differentially expressed genes (DEGs), which are involved in proteolysis, were predicted to be mediated by DEMs. These findings provide new information on the miRNA and mRNA profiles in porcine liver, which would shed light on the molecular mechanisms underlying the breed-specific traits in the pig, and may serve as a basis for further investigation into the biological functions of miRNAs in porcine liver.


Introduction
MicroRNAs (miRNAs) are short non-coding RNAs (about 22 nucleotides) which play important roles in post-transcriptional regulation by mRNA cleavage and/or translational repression [1]. miRNAs are involved in almost every biological process, including cell growth and differentiation [2], pathogenesis and disease prevention [3]. In mammals, temporal and spatial changes in miRNA expression have been well characterized [4] to delineate the miRNA transcriptomes of different tissues [2,5,6] at different development stages [7]. Also, miRNA expression profiles were elaborated on the same tissue derived from different animal breeds or species. For example, the muscle miRNA expression profiles were compared between broiler and layer chickens to understand the role of miRNAs in myogenesis [8], and miRNAs expressed in the skin of goat and sheep were profiled to study the role of miRNAs in wool growth [9].
Pig is one of the most important domestic species for meat production [10] and can also serve as an ideal model for biomedical studies on various metabolic disorders, such as obesity [11] and cardiovascular diseases [12,13] in humans. Great efforts have been made to sequence and decode the genome [14] and to identify the miRNAs [15,16] in the pig. As of November 2011, the microRNA database, miRbase (Release 18.0, http://www. mirbase.org) contains 257 mature pig miRNAs. The majority of these miRNAs are conserved among mammals, only few are pig specific. Recently, Solexa deep sequencing technology was employed not only to reveal the porcine miRNA transcriptome (microRNAome) in multiple tissues [6,7], but also to investigate the ontogeny of miRNA expression in the pig at different developmental stages [16,17].
Erhualian (EHL) is a Chinese indigenous pig breed, being known for its early sexual maturity, large litter size, high adiposity, mild temper, good maternity and high tolerance to roughage and stress [18]. In our previous study, we found distinct behavioral, endocrine and biochemical response patterns during transportation between Chinese indigenous breed and Western breed [19]. So far, mRNA transcriptomes in placenta and ovary of Chinese indigenous pigs have been profiled to understand the molecular mechanisms involved in their high prolificacy [20,21]. However, the study of hepatic mRNA transcriptome in EHL is still lacking, and the link between mRNA transcriptome and the metabolic parameters in the pig has not been well established.
Liver is a key metabolic organ and thus often being chosen as a target for miRNA profiling to understand metabolic regulations. For example, obese diabetic model mice such as ob/ob, db/db and KKAy were reported to express higher miR-335 in the liver compared to normal mice [22]. Inhibition of miR-122 in liver resulted in reduced plasma cholesterol levels, accompanied by a decrease in hepatic fatty-acid and cholesterol synthesis rate [23]. For genome-wide analysis, integrating differentially expressed miRNAs (DEMs) with differentially expressed genes (DEGs) should present a comprehensive way to study their functions in metabolomes. Some integrated analyses of liver miRNA and mRNA have been carried out in model animals [24]. However, coordinated analysis of hepatic miRNA and mRNA expression profiles in relation with the metabolic characteristics of different breeds of pigs is lacking.
Here we used male newborn piglets of Large White (LW) and EHL to delineate the expression profiles of mRNA with microarray and miRNA with both deep sequencing and microarray, and to analyze the functional relevance of integrated miRNA and mRNA expression in relation with the physiological and biochemical characters. Solexa deep sequencing was employed to get a full scope of porcine liver miRNAome and to identify the differentially expressed miRNAs between the two pig breeds. A miRNA microarray containing 238 probes and RT-qPCR were used to supplement and confirm DEMs. The mRNA expression profile was evaluated by microarray analysis. The miRNAome and transcriptome were integrated and the functional relevance was analyzed linking breed-specific metabolic phenotypes in EHL and LW pigs.

Metabolic Characteristics of Two Pig Breeds
As shown in Table 1, body weight, body length and liver weight of newborn EHL piglets were significantly lower than those of LW (FDR ,0.05), yet the liver index remained unchanged (FDR .0.05). Serum total cholesterol (TCH), high density lipoprotein cholesterol (HDLC) and low density lipoprotein cholesterol (LDLC) in EHL piglets were significantly higher than those in LW piglets (FDR ,0.05), while serum glucose and triglyceride (TG) did not differ between the two pig breeds. The TCH content in liver was also higher in EHL, while liver content of TG showed no difference between breeds. Serum cortisol concentration was significantly increased (FDR ,0.05), while serum leptin and insulin levels were significantly decreased (FDR ,0.05) in EHL piglets compared to those in LW. Serum triiodothyronine (T3), thyroxine (T4), free T3 (FT3) and T3/T4 ratio showed no difference between the two breeds (Table 1).

Liver miRNA Expression Profiling of Two Pig Breeds by Deep Sequencing, Microarray and RT-PCR Validation
After trimming of adaptor sequences and removal of reads containing ambiguous base calls, the sequence reads were clustered into unique sequences. In total, there were 25,957,969 and 26,574,154 trimmed reads in two libraries, with 15,154,209 (58.3% of trimmed reads) and 19,099,533 (71.8% of trimmed reads) mappable reads that aligned to unique miRNAs for the pooled liver samples of LW and EHL in the range of 19-25 nt. The read size mainly ranged from 21 to 23 nt. The percentage of the 22 nt reads in total reads was 64.2% for LW and 64.0% for EHL. The top 10 highest expressed miRNAs detected by deep sequencing were miR-148a, miR-101, miR-143-3p, miR-122, miR-30a-5p, miR-21, miR-30c, miR-192, miR-27b and miR-24 (Table S1).
Counts in reads per million (RPM) was used to quantify miRNA expression in the liver of EHL and LW piglets. These unique sequences with RPM .10 were annotated based on their similarities with mature miRNA sequences published in miRBase (release 18.0), resulting in a list of 211 mature miRNAs (Table S1). miRNAs with less than 2 nucleotide substitution outside seed region were considered as one miRNA family [25,26], and the copy numbers of all miRNAs in this family were added together for reporting read counts and for differential expression calculations. miRNAs with fold change .1.5 and at least 10 RPM average expression (in both pig breeds) were selected as differentially expressed miRNAs (DEMs).

Liver mRNA Expression Profiling of Two Pig Breeds by Microarray
Six pig transcriptome microarrays were performed for hepatic gene expression profiling in EHL and LW piglets. Among 44,987 pig transcripts investigated, 23,807 transcripts were retained for further analysis after removing the probes with poor quality intensities and low dependability. All these qualified transcripts were then annotated by manual blast referring to previously developed protein-based annotation for pigs [27,28]. In total, 9,447 genes were found to be expressed in the liver of newborn piglets.
The microarray data of LW piglets were treated as control in the selection of differentially expressed genes related to EHL piglets. After the removal of redundant and unannotated sequences, with FDR ,5%, 53 genes were found to be significantly up-regulated and 200 genes to be significantly down-regulated in the EHL piglets compared to the LW piglets (Fig. 1A). With FDR ,10%, 298 genes were found to be significantly up-regulated and 510 genes to be significantly downregulated in the EHL piglets (Fig. 1B).

Coordinated Expression Patterns between miRNAs and mRNAs
As described above, the DEGs were selected based on either FDR ,5% or FDR ,10% (Figure 1), while the DEMs were chosen based on either sequencing (fold change .1.5) or microarray (FDR ,15%) data.
In this study, the target gene lists of miRNAs predicted by miRanda (http://www.microrna.org/microrna/) based on the human sequences were used. Unfortunately, 10 miRNAs including 4 up-regulated (miR-129*, miR-1343, miR-193b* and miR-590) and 6 down-regulated (miR-2887, miR-2904, miR-2904-3p, miR-4332-3p, miR-4332-5p and miR-739) do not have orthologous genes in human, so they do not have target genes in the database. As such, these miRNAs were excluded for further analysis. The remaining 12 down-regulated and 26 up-regulated miRNAs were associated with 170 and 417 gene targets, respectively, in the DEG list selected at FDR ,5% (Figure 2A). Among 170 genes targeted by 12 down-regulated miRNAs, 32 (19%) were up-regulated and 138 (81%) were down-regulated. Among 417 genes targeted by 26 up-regulated miRNAs, 336 (81%) were down-regulated, and 81 (19%) were up-regulated. Similar patterns of miRNA-mRNA association were seen with the DEG list selected at FDR ,10% ( Figure 2B). Among 478 genes targeted by down-regulated miRNAs, 151 (32%) were up-regulated and 327 (68%) were down-regulated. Among 1,277 genes targeted by up-regulated miRNAs, 377 (30%) were up-regulated and 900 (70%) were down-regulated. A two tailed chi-square test was conducted to determine whether the number of predicted targets of DEMs was significantly higher than that would be expected by chance. The significance of all the 4 types of possible miRNA-mRNA correlations was analyzed. It turned out that only the correlation of up-regulated miRNAs and down-regulated genes was statistically significant (p,0.01).
To identify whether a miRNA was negatively correlated with predicted target DEGs, the number of down-and up-regulated targets of a DEM was compared with the number of stay-still targets (with FDR .10%) by two-tailed Fisher's Exact Test. The results of this analysis were presented in Figure 3. When using the DEG list of FDR ,5%, only miR-210 had a significant reciprocal expression with its targets. When using the DEG list of FDR ,10%, 10 of 26 up-regulated miRNA had a significant higher number of target mRNAs that were down-regulated. The differences obtained by using different DEG lists imply that the target of a miRNA may usually be fine-tuned. When testing miRNAs that were down-regulated, no miRNAs had a significant number of target genes that were up-regulated, but 4 miRNAs (miR-222, miR-27a, miR-574-5p, and miR-485-3p) had a significant number of target genes that were also down-regulated when using both DEG lists of ,5% and ,10%.

Functional Analysis of DEGs
To define the biological functions of all the 808 DEGs (Table  S2) selected at FDR ,10%, the gene ontology (GO) analysis and KEGG pathway analysis (Table S3) were carried out. GO terms were sorted by p-value, in an ascending order from bottom to top (Figure 4). For the biological process, liver functions between LW and EHL mainly differ in regulation of transcription, cell cycle (the GO term of cell cycle, cell division and mitosis), transcription, signal transduction, oxidation reduction, cell adhesion, lipid metabolism, DNA replication, translational elongation, development, protein amino acid phosphorylation, interspecies interaction between organisms, ion transport, immune response, response to virus, response to DNA damage stimulus, proteolysis and chromatin modification.
To characterize the function of miRNA-mediated DEGs, we used the target DEGs which are reversely expressed with DEM for GO analysis. Using the DEGs with FDR ,10%, we observed 229 coherent targets and 579 non-coherent target genes. The percentages of coherent target and non-coherent target genes involved in each GO term were and compared ( Figure 5). When considering biological process, the genes involved in proteolysis were significantly enriched among the target genes. The enrichment of proteolysis is concordant with the enrichment of hydrolase activity, peptidase activity and protein homodimerization activity in molecular function analysis. Meanwhile, DEGs involved in DNA replication, translational elongation and the constituent of ribosome are barely miRNA-mediated.

Discussion
In this study we observed significant differences in physiological and biochemical traits between newborn male EHL and LW piglets. The birth weight and liver weight indicated disparate embryo growth and organ development in the two pig breeds, as the consequence of cell proliferation and differentiation. The present results agree with our previous findings that the serum cortisol in EHL was higher than that in lean type pigs [19]. The elevated serum TCH, HDLC, LDLC and higher liver content of TCH in EHL piglets indicated higher rate of cholesterol metabolism in this breed, which was in agreement with the enrichment of lipid metabolism related genes ( Figure 4A). KEGG pathway analysis (Table S3) showed that there were 8 DEGs involved in peroxisome proliferator-activated receptor (PPAR) signal pathway. Among them, the higher expression of thrombospondin receptor (CD36) and fatty acid transporter (solute carrier family 27 member 2, SLC27A2) in EHL liver indicated increased response to serum LDLC and transport of HDLC. The lower level of CYP7A1 implicated lower synthesis of bile acid from cholesterol in the liver of EHL piglets, which may contribute to higher cholesterol deposition in liver [29]. By comparing the transcriptome of EHL and LW liver, we showed that a full scope of the liver metabolism processes varied between the two breeds of pigs.
The microarray used in our experiment could detect the miRNAs with RPM .10. For example, miR-193a-5p, which has an average RPM of 17, was detected readily by microarray. Nevertheless, the porcine miRNAs detected with deep sequencing often differ in sequence from their reference sequences in miRbase, as it was reported previously [6]. This may explain why some DEMs determined by sequencing could not be detected by microarray. The microarray probes were designed according to the miRbase sequences, whereas those abundantly expressed miRNAs (with an average RPM .10) that were undetectable by microarray had at least one nucleotide substitution compared to the miRbase sequences. Occasionally, some miRNAs with mismatches could also be detected by microarray. For example, miR-100, which has an average RPM of 958, was detected, however with the hybridization signal much lower than the perfectly matching miRNAs of the similar abundance. Besides, low copy miRNAs might also be undetectable by microarray due to low hybridization signals. Therefore, we combined the informa-tion obtained from deep sequencing and microarray formats to maximize the list of DEMs between EHL and LW piglets.
The coherent relationship between miRNA and mRNA is still under debate. Initial study showed that some miRNAs could induce reversed expression of mRNA and protein [30]. Actually, some recent studies showed that the uncoupling between mRNA and protein may implicate the post-transcription regulation mechanism [31]. Till now, there are 110 examples of mRNA cleavage and 178 examples of mRNA translation repression in Tarbase 5.0 (experimentally proved miRNA targets database) [32,33]. But all the high-throughput experimental methods for identifying miRNA targets usually identify mRNAs or proteins which are down-regulated when a miRNA is over-expressed or vice versa [34]. It was indicated that the miRNA-induced destabilization of target mRNAs is the predominant reason for reduced protein output [35]. In agreement with the previous studies, we found that only the target genes with decreased expression corresponding to the up-regulated miRNAs were significantly enriched. Although 4 up-regulated miRNAs were found to be associated with the up-regulated targets (p,0.05), the correlation of up-regulated miRNAs and up-regulated mRNA was not statistically significant in general. So we integrated the reversely expressed targets of DEMs to delineate the regulatory mechanisms of miRNAs between EHL and LW piglets.
Most miRNAs in mammals pair imperfectly with their target mRNAs. Therefore it is difficult to seek their biologically important targets by the algorithm analysis alone [36]. Many algorithms are based on seed pairing-the paring of miRNA nucleotides 2-8. The miRanda, which was developed by miRbase, had the latest information of miRNA targets. Owing to the insufficiency of pig gene sequence annotation, the untranslated regions (UTRs) of many pig genes have not been identified. That is why we used the human genome information to generate our miRNA target list in the present study.
Compared to LW piglets, EHL piglets had more up-regulated miRNAs than down-regulated miRNAs in liver, as revealed by both deep sequencing (25 up-regulated, 15 down-regulated) and microarray (7 up-regulated, 3 down-regulated) methods. At the mRNA expression level, the figure went to opposite, i.e., there were less upregulated than down-regulated mRNA genes (53 up-regulated, 200 down-regulated with FDR ,5%, and 298 up-regulated, 510 downregulated with FDR ,10%). miRNAs are known to exert posttranscriptional regulation mostly by inducing mRNA degradation. Therefore, miRNAs are often expressed in an opposite pattern to the mRNA expression level of their target genes. Such inverse correlations between the expression of miRNAs and their target mRNAs are also observed in other studies [37,38].
In an earlier study, the miRNA expression in liver was compared between Tongcheng (another Chinese indigenous breed) and Large White pigs at about 25 kg body weight by microarray [6]. Forty five miRNAs were found to be up-regulated, while only 13 were downregulated in Tongcheng pigs. Among these DEMs, miR-133a, miR-451 and miR-739 are also identified as DEMs in the present study. The predominant up-regulated miRNA expression in the liver of Chinese indigenous pig breeds represented that the genesis of miRNA in these breeds was more active than that in LW. In fact, the enzymes involved in microRNA processing, Drosha and Dicer, were significantly up-regulated in EHL piglets than that in LW at protein level (data not shown). The miRNA mediated processes are identified by analyzing the function of the target DEGs which were reversed correlated with DEM. Majority of the genes involved in DNA replication and protein translation were predicted to be the non-target genes of miRNA regulation (Figure 5), possibly because that genes participating in the DNA duplicating process or consisting ribosomes (ribosome proteins and histones) usually have short 39UTR [39]. In the present study, the role of miRNAs in mediating the process of lipid metabolism was not evident, despite the fact that the differences in lipid metabolism are significant between the two pig breed. Nevertheless, miR-335, which was previously shown to be expressed more abundantly in the liver of obese mice [22], demonstrated higher hepatic expression in EHL piglets. In contrast, proteolysis is predicted to be a miRNAmediated process. There are 28 DEMs targeting all the 14 genes involved in proteolysis, indicating that the process of proteolysis is regulated by the cooperation of many miRNAs. Among the 14 DEGs, 11 were down-regulated while 3 were up-regulated in EHL piglets. The decreased expression of most proteolysis related genes (11 of 14) in the liver of EHL indicated a weaker ability of protein turnover, and probably lower growth rate. However, mRNA expression is not directly related to the function. In some cases, the levels of mRNA and protein are reversely correlated. Therefore, it is possible that some proteolysis related DEGs are upregulated (3 of 14) while the majority of the relevant genes are down-regulated. Nevertheless, direct evidences regarding the breed differences in hepatic proteolysis are needed to support our presumption.
In conclusion, we demonstrated the differences in the hepatic transcriptome and miRNAome between EHL and LW piglets with distinct phenotype and metabolic character. The most highly miRNA-mediated biological process with significant breed disparity is proteolysis. Our findings provide new information on the miRNA and mRNA profiles in porcine liver, and new approach in characterizing diverse traits in different pig breeds, thus serving as a basis for further investigation of the biological functions of miRNAs in porcine liver.

Animal Sampling
The newborn piglets were obtained from two neighboring pig breeding farms and sacrificed immediately after birth by exsanguination. The experiment was conducted following the guidelines of Animal Ethics Committee at Nanjing Agricultural University, China. The slaughter and sampling procedures complied with the ''Guidelines on Ethical Treatment of Experimental Animals' ' (2006) No. 398 set by the Ministry of Science and Technology, China and ''the Regulation regarding the Management and Treatment of Experimental Animals'' (2008) No.45 set by the Jiangsu Provincial People's Government. Six newborn male piglets from three litters (2 from each litter) of each purebred EHL and LW sows were sacrificed. Body weight, body length and liver weight were recorded. The blood was collected from the precaval vein and the serum was gathered and kept at 220uC. Liver samples were immersed in liquid nitrogen immediately after collection and then stored at 270uC.

Measurement of Serum Biochemical Parameters
The measurement of serum glucose, TG, TCH, HDLC and LDLC was carried out by a service provider, Nanjing Jiancheng Bioengineering Institute.
All data were expressed as mean 6 SEM. The p-value was calculated by the student t-test for independent samples with the SPSS 13.0 for Windows. False discovery rates (FDRs) were calculated using the method of Benjamini & Hochberg (1995) [40]. FDR ,0.05 was considered significant.

RNA Isolation
Total RNA was isolated from liver,using the Trizol reagent (Invitrogen, USA), according to the manufacturer's instructions. Concentration of the extracted RNA was measured using NanoDrop ND-1000 Spectrophotometer [41]. RNA integrity was confirmed by denaturing agarose electrophoresis, and DNA contamination was evaluated by PCR using isolated RNA as template with the primers of GAPDH (primer sequences are shown in Table S4). High quality RNA samples were then stored at 270uC till further use.

Small RNA Library Preparation and Sequencing
The total RNA samples were prepared as follows: 2000 ng of total RNAs from the twelve pigs were isolated and pooled within each breed (EHL and LW). Both small RNA libraries were generated according to Illumina's sample preparation instruction. Then they were sequenced on the Illumina GAIIx (Illumina, USA) following vendor's instruction. Raw sequencing reads were obtained using Illumina's Pipeline v1.5 software following sequencing image analysis by Pipeline Firecrest Module and base-calling by Pipeline Bustard Module.
The reads were then subjected to a series of data filtration steps to obtain mappable sequences using ACGT101-miR v3.5 (LC Sciences, USA) with the statistics of mammalian miRNAs in miRBbase 16.0.

miRNA Microarray
Total RNAs of the two littermate piglets were pooled for the microarray analysis. The pooled samples were named as EHL1-3 (n = 3) and LW1-3 (n = 3).
The pig microRNA microarray was obtained from LC Sciences (Houston, USA) and contains 238 unique probes that were complementary to all mature miRNAs of pig in miRBase release 16.0. The assay started from 5 mg total RNA sample, which was size fractionated using a YM-100 Microcon centrifugal filter (Millipore, USA) and the small RNAs (,300 nt) isolated were 39extended with a poly(A) tail using poly(A) polymerase. An oligonucleotide tag was then ligated to the poly(A) tail for later fluorescent dye staining. Hybridization was performed overnight on a mParaflo microfluidic chip using a micro-circulation pump (Atactic Technologies, USA). On the microfluidic chip, each detection probe consisted of a chemically modified nucleotide coding segment complementary to target miRNAs (from miRBase, release 16, 238 probe sets) and a spacer segment of polyethylene glycol to extend the coding segment away from the substrate. The detection probes were made by in situ synthesis using photogenerated reagent (PGR) chemistry. The hybridization melting temperatures were balanced by chemical modifications of the detection probes. 100 mL 66SSPE buffer containing 25% formamide was used for hybridization at 34uC. After hybridization, fluorescence labeling using tag-specific Cy5 dyes was used for detection. Hybridization images were collected using a laser scanner (GenePix 4000B, Molecular Device) and digitized using Array-Pro image analysis software (Media Cybernetics). Data were analyzed by first subtracting the background and then normalizing the signals using a LOWESS (locally weighted regression) function. The three samples of LW piglets were used as control group to analysis the different expression between the two pig lines. The differences between the two groups were analyzed using SAM (Significance Analysis of Microarrays) method with SAMR software version 3.02 [42]. The FDR (False Discovery Rate) and fold change were calculated. miRNAs with FDR ,15% were considered to be differentially expressed miRNAs. All the microarray data were MIAME compliant and have been deposited in GEO (accession number GSE33523).

Real-time qPCR Confirmation of Different Expressed miRNAs
Two mg of total RNA from each piglet were polyadenylated by poly(A) polymerase (PAP) at 37uC for 1 h in a 20 mL reaction mixture using Poly(A) Tailing Kit (AM1350, Ambion, USA) according to the manufacturer's instructions. Treated RNA was then dissolved and reverse-transcribed using poly (T) adapter and M-MLV (Promega, USA). qPCR for the 18 miRNAs was performed using SYBR Green Real-time PCR Master Mix (TaKaRa, Japan) in Mx3000P (Stratagene, USA) with a miRNA specific forward primer and a universal reverse primer complementary to part of the poly(T) adapter sequence. U6 was chosen as an internal control to normalize the technical variations. The sequences for all the primers and the poly (T) adapter are listed in Table S4. The method of 2 2DD CT was used to analyze the realtime PCR data expressed as the fold change relative to the average value of the LW piglets [43]. The differences between the two groups were analyzed by the student t-test for independent samples with the SPSS 13.0 for Windows.

mRNA Microarray Experiment
The samples used for mRNA microarray were the same as the miRNA microarray. Microarray was performed by a service provider (CapitalBio, China). The microarray (Probe length 60mer, 135K Format) containing 44987 probe sets was provided by Roche-NimbleGen. For each sample, 1 mg of total RNA was treated with the CapitalBio cRNA Amplification and Labeling Kit (CapitalBio, China) according to the manufacturer's instructions. After reverse transcribed with a T7 oligo(dT) primer, secondstrand synthesis and purification, the generated cDNAs were in vitro transcribed to synthesize multiple copies of cRNAs. Then 5 mg of purified cRNAs were reverse transcribed with random primer. Labeled cDNA molecules were generated by subsequent Klenow Fragment Polymerase labeling with Cy3-dCTP (GE Healthcare, USA). Following purification and drying, the labeled cDNAs were then dissolved in 80 ml hybridization buffer (36SSC, 0.2%SDS, 56Denhart's, 25% formamide). Hybridizations were performed overnight at 42uC using hybridization system 12 (Roche NimbleGen, USA). The arrays were then washed and dried. The fluorescence intensity was collected using NimbleGen MS 200 Microarray Scanner. Data were extracted from scanned images using NimbleScan v2.6 software. Quantile normalization RMA (Robust Multi-Array) analysis was performed to generate gene expression values. The differences between the two groups were analyzed using SAM (Significance Analysis of Microarrays) method with SAMR software version 3.02 [42]. The FDR (False Discovery Rate) were calculated. Differentially expressed genes (DEGs) were selected with FDR ,5% and FDR ,10%. All data were MIAME compliant and have been deposited in GEO (accession number GSE33524).

Bioinformatics and Statistical Analysis
The target genes of miRNAs were predicted by miRanda algorithm. Correlation analysis of the miRNA and mRNA expression profiles was carried out using the lists of DEMs and DEGs. The significance of all the 4 types of possible miRNA-mRNA correlations (up-regulated miRNA and up-regulated mRNA, up-regulated miRNA and down-regulated mRNA, down-regulated miRNA and down-regulated mRNA, downregulated miRNA and down-regulated mRNA) were analyzed using two tailed chi-square test.
A two-tailed Fisher Exact Test was conducted for each DEM to determine whether the number of predicted target genes that were differentially regulated was higher than would be expected by chance. The Fisher Exact test was conducted for each miRNA using both down-regulated and up-regulated DEG lists.
The GO and KEGG pathway analysis were carried out by using a Molecule Annotation System called MAS 3.0 (http:// bioinfo.capitalbio.com/mas3/) and the enrichment p-values were calculated.
The DEG list with FDR ,10% was used to characterize the function of miRNA-mediated DEGs. We classified the DEGs into two categories, coherent targets and non-coherent targets. Coherent target genes are predicted DEM target genes that are negatively correlated with the expression of DEMs. For each GO term, the percentages of the coherent and non-coherent targets were compared using a two-tailed chi-square test.