miRNA Expression Profile Analysis in Kidney of Different Porcine Breeds

microRNAs (miRNAs) are important post-transcriptional regulators in eukaryotes that target mRNAs repressing their expression. The uncertain process of pig domestication, with different origin focuses, and the selection process that commercial breeds suffered, have generated a wide spectrum of breeds with clear genetic and phenotypic variability. The aim of this work was to define the miRNAs expression profile in kidney of several porcine breeds. Small RNA libraries from kidney were elaborated and high-throughput sequenced with the 454 Genome Sequencer FLX (Roche). Pigs used were classified into three groups: the European origin group (Iberian breed and European Wild Boar ancestor), European commercial breeds (Landrace, Large White and Piétrain breeds) and breeds with Asian origin (Meishan and Vietnamese breeds). A total of 229 miRNAs were described in the pig kidney miRNA profile, including 110 miRNAs out of the 257 previously described pig miRNAs and 119 orthologous miRNAs. The most expressed miRNAs in pig kidney microRNAome were Hsa-miR-200b-3p, Ssc-miR-125b and Ssc-miR-23b. Moreover, 5 novel porcine miRNAs and 3 orthologous miRNAs could be validated through RT-qPCR. miRNA sequence variation was determined in 116 miRNAs, evidencing the presence of isomiRs. 125 miRNAs were differentially expressed between breed groups. The identification of breed-specific miRNAs, which could be potentially associated to certain phenotypes, is becoming a new tool for the study of the genetic variability underlying complex traits and furthermore, it adds a new layer of complexity to the interesting process of pig evolution.


Introduction
MicroRNAs (miRNAs) are small regulatory RNAs that play important roles in the regulation of gene expression [1,2]. These single-stranded non-coding RNAs that are approximately 22 nucleotides long are involved in post-transcriptional regulation mechanisms acting mainly through down-regulation of target messenger RNAs (mRNAs) in a wide range of biological and pathological processes [3,4]. Thus, miRNAs inhibit gene expression by blocking protein translation or inducing mRNA degradation [5][6][7]. mRNAs can be regulated by several miRNAs, and each miRNA can target hundreds of mRNAs in different binding sites [8,9]. The recent evidence about the genetic variability at both 59 and 39 ends of the mature miRNA sequence generates a large spectrum of sequence variants called miRNA isoforms or isomiRs [10][11][12]. Mechanisms involved in isomiRs generation and their biological relevance has increased the complexity of molecular mechanisms related to regulation of mRNA expression mediated by miRNAs.
Nowadays, a large number of miRNAs have been reported in animals, plants and viruses, with up to 18,226 entries in the miRBase database (v18, November 2011, URL: http://www. mirbase.org/) [13][14][15]. In mammalian genomes, the number of encoded miRNAs has been predicted to be up to 1000 miRNAs, comprising about 3% of all protein-coding genes [5]. In mammals, nearly all miRNAs are conserved in closely related species [16], and may have homologous in distant species, suggesting that miRNA functions could also be conserved throughout the evolution of animal lineages [17]. Several studies showed that variability in miRNA sequences has been lost over time since miRNAs have been described in unicellular eukaryotes, showing its deeper evolutionary history among eukaryotes [18].
The pig is an important livestock species, not only for its production in industry, but also as a suitable animal model for comparative genomics and biomedical studies [19]. However, the number of porcine miRNAs available in public databases is still poor, with only 257 porcine miRNAs described in the pig genome compared to the completely sequenced human (1,921), mice (1,157), bovine (676), poultry (544), equine (360) or canine (289) genomes (miRBase v18). The description of the microRNAome (miRNAome) held by different tissues under different pathophysiological states has become of great interest in the last years [20][21][22][23][24]. Determining altered patterns of miRNA expression related to disease or specific treatments or conditions would be useful in order to identify differentially expressed miRNAs that could be used as novel biomarkers [25][26][27][28][29][30]. In pigs, several miRNAomes have been described in different tissues such as muscle, fat, heart, liver, thymus, intestine and testes [31][32][33][34][35], or using in vitro cells models such as porcine PK-15 cells (derived from porcine kidney epithelial cells) and porcine dendritic cells [36][37][38]. However, the porcine kidney miRNAome has not yet been described although it is an essential organ involved in functions like blood filtering, gluconeogenesis and in the secretion of important hormones like erythropoietin, renine and vitamin D.
The process of pig domestication has been very complex. Independent geographical origins have been described [39] which have generated multiple phenotypically different breeds [40,41]. In particular, high differences are found regarding reproductive and meat production and quality traits between Asian and European breeds. In addition, pig has become an important production animal for human, and, consequently, it has been strongly selected for traits of economical interest. As a result of this high selection, many commercial breeds have been generated. Different consuming necessities in pig industry have expanded the large phenotypic variability due to different selected traits. Genetic diversity between pig breeds has been further studied [42][43][44] showing variability in gene expression and elucidating the artificial selection performed by porcine industry during pig domestication [45]. However, the genetic variability held at miRNA sequences, which could imply changes at post-transcriptional level, has not yet been explored although it may lead to great phenotypic differences and/or pathological disorders.
Studies with mitochondrial DNA (mtDNA) assume significant differences between European and Asian pigs prior to domestication and several studies describe the phylogenetic relationships between European and Asian pig breeds, determining different groups into European diversity and other groups into Asian diversity [39,44,46]. Other studies confirmed some introgression events of Asian pig breeds in most European commercial breeds and evidence that local Spanish breeds have never been introgressed with Asian neither with European commercial populations [47][48][49]. In this sense, studying miRNAs differentially expressed between European and Asian pig breeds could contribute to understand their phenotypical differences.
In the present work, the miRNAome of several pig breeds belonging to the pure European and Asian branches as well as European commercial breeds that have been introgressed with Asian genes has been undertaken. The main goals of this study were to describe the pig kidney miRNAome through high throughput sequencing (HTS) technology, to study the miRNA sequence variability among porcine breeds and to determinate new porcine miRNAs.

Samples and RNA preparation
Kidney samples were collected from twenty pigs belonging to six different breeds and the European ancestor (Wild boar). Animals were classified into three main groups based on breed origin; (1) European breeds: Wild boar (WB, n = 2) and Iberian (IB, n = 4), (2) Asian breeds: Meishan (ME, n = 3) and Vietnamese (VT, n = 3), and (3) European commercial breeds with strong influence from Asian breeds: Landrace (LD, n = 3), Large White (LW, n = 2) and Piétrain (PT, n = 3). All samples were taken from slaughterhouse (PRIMAYOR, Mollerussa, Spain) under veterinary supervision and they were immediately snap-frozen in liquid nitrogen and stored at 280uC until use. Total RNA was isolated using TRIzolH reagent (Invitrogen, Carlsbad, USA) following the manufacturer's recommendations. Total RNA was quantified using ND 1000 NanodropH Spectrophotometer (Thermo Scientific, Wilmington, USA) and its quality was assessed on an Agilent 2100 Bioanalyzer using the RNA 6000 Nano kits (Agilent Technologies, Santa Clara, USA). RNA quality threshold was set to RNA Integrity Number (RIN) above seven. For RT-qPCR validations, additional kidney samples were added (n = 29: 2 WB, 4 IB, 4 LD, 4 LW, 4 PT, 4 ME and 7 VT).

Small RNA library construction and high throughput sequencing
Small RNA fraction (15-30 nt) was excised and isolated from denaturing 12.5% polyacrilamide gel electrophoresis (PAGE) using miSpike TM (IDTH, Coralville, USA) as internal size marker. For each breed, 60 mg of total RNA were loaded on separate gels to avoid cross-contamination. Gels were stained with GelStarH Acid Nucleic Gel Stain (Lonza, Basel, Switzerland) for UV visualization. Excised small RNA fraction were purified using PerformaH DTR gel filtration cartridges (EdgeBio, Gaithersburg, USA) and, then, isolated small RNAs were pooled by breed to construct small RNA libraries. Briefly, 39 and 59 linkers from miRCat TM kit (IDT, Coralville, USA) were ligated at both ends of the small RNAs in two separated reactions using a T4 RNA ligase without ATP (Fermentas, Germany) and T4 RNA ligase with ATP (Ambion, Austin, USA), respectively. Between 39 and 59 primer ligations, the 60 nt RNAs were purified by PAGE to eliminate unligated products. Then, linked products were used to perform a reverse transcription reaction using the SuperScript TM III Reverse Transcriptase kit (Invitrogen TM , Carlsbad, USA) and the cDNA obtained was amplified with the PlatinumH Taq DNA Polimerase High Fidelity kit (Invitrogen TM , Carlsbad, USA). PCRs were done with primers complementary to 39 and 59 linkers and, in addition, they included multiplex identifiers at the 59 end (a five nucleotide sequence tag) to allow differentiation between libraries (Table S1). Obtained PCR products were purified using the QIAquick PCR Purification Kit (QiagenH, Germany). Libraries were quantified with Qubit TM fluorometer, Quant-IT TM (Invitrogen TM , Carlsbad, USA) and prepared to a 10 8 DNA molecules/mL concentration for sequencing by 454 Genome Sequencer FLX Titanium (Roche, Germany) following the manufacturer's protocol at DNA sequencing facilities at CRAG (Bellaterra, Spain).
Sequence annotation for miRNA expression profile and validation of selected miRNAs through Reverse Transcription quantitative real time PCR (RT-qPCR) Primers were trimmed from sequences and only those sequences between 15 and 29 nucleotides and with total number of sequences $3 were kept for further analysis. Sequences were compared to all available miRNA sequences (miRBase v18) using local Blast. Parameters were set to 100% identity and up to 4 mismatches allowed at the end of the sequences to assume variability on 39 and 59 ends [31].
Reverse transcription (RT) reactions were performed in duplicate with the universal cDNA synthesis kit (Exiqon, Denmark) using 300 ng of total RNA following the manufacturer's instructions. Non template controls (NTC) and minus poly(A) polymerase controls for each tissue were included. qPCRs were performed in duplicate using miRCURY TM LNA TM Universal RT microRNA PCR kit (Exiqon, Vedbaek, Denmark) on an 7900HT Sequence Detection System (Applied Biosystems, Warrington, UK). Each qPCR was done in 20 mL total volume including 10 mL SYBR Green master mix (Exiqon, Vedbaek, Denmark), 0.5 mL of each LNA primer set (Exiqon, Denmark) and 5 ml of a 1:20 dilution of the cDNA. To assess qPCR efficiency, standard curves with 10-fold serial dilutions of a pool of equal amounts of cDNA from all samples were included in each assay. Thermal profile was set as follows: 95uC for 10 min and 40 cycles at 95uC for 15 sec and 60uC for 60 sec. NTC and minus poly(A) polymerase controls were included. Melting curve analysis was included at the end of the qPCR to detect unspecific amplifications.
Quantities from each sample were obtained from the calibration (standard) curve added in each RT-qPCR reaction. GeNorm v.3.5 software [53] was used to examine the stability of the reference miRNAs (M,1.5) and to obtain a normalization factor (NF). The quantity obtained from each miRNA and sample was normalised by the NF and FCs were calculated in relation to the lowest normalised value. FCs were log 2 transformed and expression data were analysed by a two-way analysis of variance (ANOVA) with the General Linear Models procedure of the Statistical Package for the Social Scientists (IBMH SPSSH Statistics 19; IBM Corporation, Armonk, USA) including the RT and breed as fixed factors. Significance threshold was set at a,0.05. Scheffe test was used to determine significant differential expression between breed groups.
Sequence analysis for novel miRNAs discovery and validation through RT-qPCR All primer trimmed sequences were mapped in the pig genome retrieved from Ensembl Genome Browser (Ensembl release 62, ftp://ftp.ensembl.org/pub/release-62/fasta/sus_scrofa/dna/, Sscrofa9.62) considering 100% of alignment and identity (perfect match). Sequences that mapped only once in the genome and with a length from 19 to 23 nt were selected. Among these, only sequences with unknown annotation at Ensembl genome browser and with copy number (CN) higher than 2 were considered. They were clustered taking into account only the position in the genome. Hence, sequences positioned in the same region were grouped and the sequence with higher CN was selected as the reference sequence for each cluster. Fifteen clusters were considered novel miRNAs to be validated. Flanking regions (50 nt) of the selected reference sequences were used to predict pre-miRNA folding structure using MFold software [54] following the guidelines reported by Ambros et al. [55]. At the end, 8 clusters were selected for RT-qPCR validation. RT reactions were performed in duplicate using total RNA as previously described by Balcells et al [52]. Briefly, 600 ng of total RNA in a final volume of 20 mL including 2 mL of 10x poly(A) polymerase buffer, 0.1 mM of ATP, 0.1 mM of each dNTP, 1 mM of RT-primer, 200 U of M-MuLV Reverse Transcriptase (New England Biolabs, USA) and 2 U of poly(A) polimerase (New England Biolabs, USA) was incubated at 42uC for 1 hour and at 95uC for 5 minutes for enzyme inactivation. Non template controls (NTC) and minus poly(A) polymerase controls for each tissue were included.
DNA primers for each miRNA were designed following the methodology suggested by Balcells et al [52] (Table 1). qPCR reactions were performed in duplicate in 20 mL final volume including 10 mL FastStart Universal SYBR Green Master (Roche, Germany), 125-500 nM of each primer (Table 1) and 5 mL of a 1:100 dilution of the cDNA. Standard curves were generated by 10 fold serial dilutions of a pool of all cDNAs in order to calculate the qPCR efficiency. For Cluster 2 (Cl-2), standard curve was performed by using 2-fold serial dilutions. qPCR settings and data analysis was performed as explained above. Hsa-let-7a, Hsa-miR-25 and Hsa-miR-103 were used as reference miRNAs [50,51]. Target prediction and functional analysis DIANA -microT v3.0 web server [56,57] was used to identify in silico potential mRNA targets for differentially expressed miRNAs. Porcine genes are not included in the current version of DIANA -microT v3.0 and predictions were based on the human mRNA: miRNA interactions assuming sequence conservation. In silico functional annotation of putative mRNA target genes for each miRNA were analyzed with WEB-based Gene Set Analysis Toolkit (WebGestalt, [58]). Predicted miRNA targets were functionally annotated through the biological process information supported by Gene Ontology (GO, [59]) and the pathways in which they were involved were described by using the Kyoto Encyclopedia of Genes and Genomes database (KEGG, [60,61]). Over or under represented functional categories were identified with hypergeometric test corrected by the multiple test adjustment proposed by Benjamini & Hochberg [62]. Significant threshold was set at a,0.05.

Characterisation of miRNA expression profile in porcine kidney
Small RNA libraries from kidney of six porcine breeds (Iberian, Landrace, Large White, Piétrain, Meishan and Vietnamese) and the European ancestor (Wild Boar) were sequenced in a high throughput 454 GS FLX sequencer (Roche). A total of 115,305 reads (corresponding to 17,913 unique sequences) containing short inserts of length ranging from 15 to 29 nt (corresponding to miRNA size) were obtained ( Table 2). The total number of reads obtained per breed ranged from 13,704 to 23,815 counts; however, the Iberian breed yielded a much lower number of total counts (4,022). If we consider the number of unique sequences obtained per breed, as an average, 17% of total counts corresponded to unique sequences except for Iberian and Vietnamese breeds that had 39% and 38% of unique sequences, respectively.
Alignment of sequences to miRBase database revealed that 77.7% (89,594) of the total sequences yielded a positive match to a known miRNA. In total, 229 different miRNAs were found to be expressed in porcine kidney. Among them, 110 miRNAs corresponded to previously described porcine miRNAs and 119 to orthologous miRNAs. The percentage of detected miRNAs per breed varied from 55 to 83%, being the Iberian breed the library with less detected miRNAs and the Large White the most, in agreement with the total number of sequences obtained in these libraries. The proportion of porcine miRNAs and orthologous miRNAs was similar in all libraries, around 50% both. Sixty-two miRNAs had a CN below 3 and were excluded from the analysis. Overall, the porcine kidney miRNA profile was estimated to be of 167 miRNAs (Table S2). The most expressed miRNAs (CN.350, 0.30%) in porcine kidney are listed in Table 3. The most abundant miRNA was Hsa-miR-200b-3p (27,097, representing 23.50% of all porcine kidney miRNAs), followed by Ssc-miR-125b (8,809; 7.64%), Ssc-miR-23b (5,412; 4.69%), Ssc-miR-126 (5,274; 4.57%) and Ssc-miR-23a (5,156; 4.47%). The miRNA expression pattern was similar between breeds; however, some differences were found. As an example, Ssc-miR-125b was the most expressed miRNA in Iberian breed, followed by Ssc-miR-99a and Hsa-miR-200b-3p, and Ssc-miR-192 and Hsa-miR-200c-3p were the fourth and fifth most expressed miRNAs in Large White breed, respectively. In addition, Hsa-miR-200c-3p appeared as the fifth most expressed miRNA in Meishan breed.

Sequence miRNA variation
All miRNA sequences showed variability at both ends although it was higher in the 39 end. In some cases, nucleotide variability could be also detected inside the miRNA sequence. In general terms, for each miRNA there was one sequence more expressed called reference miRNA sequence and, then, there was a number of sequence variations with less number of copies ( Figure 1). Nucleotide variation was found in 70% of the miRNAs expressed in kidney and only a small proportion of kidney miRNAs did not show isomiRs. Nevertheless, those miRNAs without sequence variations were low expressed, from 3 to 38 total reads. Overall, 116 miRNAs presented isomiRs with an average rate of 10 isomiRs per miRNA. The miRNA with more variants was Hsa-miR-200b-3p with 123 variants, followed by Ssc-miR-23b and Ssc-miR-125b, with 59 and 51 variants, respectively. These three miRNAs are, in fact, the most expressed miRNAs in porcine kidney (Table 3). In this sense, variants expression followed two main patterns: according to those miRNAs with more than 1,000 total reads, they were distributed in those miRNAs with a strong predominant isomiR, such as Hsa-miR-200b-3p, Ssc-miR-125b, Ssc-miR-23b, Ssc-miR-23a, Ssc-miR-192, Ssc-miR-10b, Ssc-miR-126* and Ssc-miR-10a, and those miRNAs where there is not a really strong predominant isomiR, like Ssc-miR-126, Ssc-miR-99a, Hsa-miR-200c-3p, Ssc-miR-30d and Ssc-miR-125a (Table S3). Furthermore, the most expressed isoform was, in most cases, the same in all breeds; it only differed between breeds in Ssc-miR-99a and Has-miR-200c-3p.
Interestingly, reference miRNA sequence for each miRNA was not always the same as the one described in miRBase database. This was the case of 72 miRNAs (Table S4). Sequence differences were mainly found at 39 end of the miRNA (in 90% of cases) and consisted, basically, in the addition or deletion of one or two nucleotides. In Table 4 are shown those most expressed miRNA sequences (CN.100) differing from the miRBase reference sequence. These differences in miRNA sequences are not influenced by the miRNA abundance, as total reads in these 72 miRNAs include high expressed miRNAs (Ssc-miR-23a, 5,156 reads) and low expressed miRNAs (Mmu-miR-5100, 3 reads). Moreover, counts of these 72 miRNAs represent the 21% of total counts of kidney miRNA profile, a considerable amount.
Furthermore, the most abundant reference sequence for a specific miRNA varied also between breeds. We found 76 miRNAs for which the most abundant sequence in each breed was not the same (data not shown). Nonetheless, differences in isomiR predominance between breeds did not follow any pattern related to breed origin. In addition, for some miRNAs, the described sequence in miRBase was expressed at a low level or even could not be detected in the porcine kidney samples analysed in this study. This was the case of Bta-miR-193b, Ssc-miR-374a, Ssc-miR-145, Ssc-miR-362, and Hsa-miR-324-3p, that had a CN of 21, 13, 3, 0 and 0, respectively.
To propose a nomenclature for the novel miRNAs, they were aligned to miRBase database accepting internal mismatches and some variation at the ends (Table 6). Cl-5 had a positive match with Hsa-miR-194-5p (p-value = 0.0004) with two nucleotide deletions at the ends, and also with Ssc-miR-194 (p-value = 0.002) although in this alignment one internal mismatch was detected in addition of the two nucleotide deletions at the ends. Cl-25 had homology with Hsa-miR-551a and Cl-29 did not match to any orthologous miRNA in miRBase database. Cl-38 matched with Bta-miR-1468 (p-value = 0.0003) with one nucleotide deletion in the 59 end. Cl-2 was not differentially expressed through RT-qPCR but it matched with Hsa-miR-31-3p (p-value = 0.001) having two nucleotide deletion at the ends.

Target prediction and functional analysis of differentially expressed porcine miRNAs
Target genes for the eight differentially expressed miRNAs (Hsa-miR-200b-3p, Hsa-miR-200c-3p, Ssc-miR-126, Ssc-miR-126 * , Ssc-miR-99a, Bta-miR-193b, Ssc-miR-486 and Ssc-let-7f) were predicted in silico. A total of 1,884 target genes were identified (Table S8) and they were functionally analyzed through KEGG pathways database where 96 biological pathways were determined (data not shown). Interestingly, Hsa-miR-200b-3p and Hsa-miR-200c-3p, from the same miRNA family, practically share the same pathways, related to cellular processes, signal transduction and some biological processes, like immune, nervous, circulatory and endocrine systems. They have also been related to cancer pathways, like renal cell carcinoma, among others. Ssc-miR-126, Ssc-miR-126*, Ssc-miR-99a and Ssc-let-7f have also been described in some cancer pathways, emphasizing the importance of miRNAs in pathways related with cell cycle, growth and death. No significant related pathways were obtained from targets associated to Bta-miR-193b and Ssc-miR-486.

Discussion
The present study represents the first work that describes the porcine kidney microRNAome in 6 different breeds and the European Ancestor, the Wild Boar, to determine differentially expressed miRNAs related to breed origin. A total of 110 out of 257 known porcine miRNAs in miRBase database have been determined in this study by high throughput sequencing. It is known that some miRNAs could be tissue specific [64] and this could explain why this study has only described the 43,28% of the porcine miRNAs reported up to date which is in agreement with other studies [22,65]. Interestingly, among all described kidney miRNAs, 119 mammalian orthologous miRNAs were identified. Assuming high miRNA conservation in mammals [16,66], the Figure 1. Sequence variability in Ssc-miR-23b. The number of total counts for each sequence is indicated. First sequence belongs to the precursor miRNA. Red bases match with precursor miRNA and belong to mature miRNA. Blue bases match also with precursor miRNA but do not belong to mature miRNA sequence. Sequence corresponding to the annotated miRNA in miRBase database is marked with an asterisk. doi:10.1371/journal.pone.0055402.g001 high proportion of orthologous miRNAs (52%) suggests that there are still many miRNAs not described in pigs yet. miRNA expression profile in kidney revealed that the most expressed miRNAs were Hsa-miR-200b-3p, Ssc-miR-125b and Ssc-miR-23b. Interestingly, Hsa-miR-200b-3p is the most expressed miRNA in all breeds except for Iberian breed, which Ssc-miR-125b (567 reads), Ssc-miR-99a (339 reads) and Hsa-miR-200b-3p (239 reads) are the first, the second and the third most expressed miRNAs, respectively. This variation of the expression in Iberian breed could be explained by many reasons. Firstly, some bias could be assumed because of the low number of reads obtained in this library in comparison with the other libraries (Table 2). However, the total reads of each miRNA in Iberian breed evidence that there is some tendency in its expression. Secondly, this change in the most expressed miRNA could also be explained as Iberian breed was exclusively originated from European Wild Boar and not mixed with any other ancestor neither being in contact with other commercial breeds in Europe.
Importantly, the most expressed miRNA in kidney was Hsa-miR-200b-3p which has not already been described in pig. Its sequence could not be found in the previous pig genome sequence (v9.2) but it could be now annotated in the last version (Sscrofa10, ftp://ftp.ensembl.org/pub/release-67/fasta/sus_scrofa/dna/) probably due to it was located on a gap in the pig genome sequence that has been solved in the new version. Its expression was significantly higher than other profiled miRNAs, as it had been detected 27,097 times corresponding to the 23.50% of the total reads. This expression pattern in miRNAs profile has been defined in more studies [67], where a miRNA is predominant over the rest meaning that it could have a significant role in the expressed tissue. It is hypothesised that it could governs or be implicated on the major constitutive functions carried out by this tissue.
The multiple sequences obtained for each specific miRNA were aligned to study the variability in the 39 and 59 miRNA ends. Most of the identified miRNAs (70%) showed nucleotide polymorphisms and variations in sequence length, according to previous studies [31,68]. One example in this study is Ssc-miR-23b, where the length varies from 18 to 25 nucleotides (Figure 1). This variability in the 59 and 39 ends could be favoured by processes prior to the high throughput sequencing sample preparation steps [69], being explained by the cleavage of the same miRNA precursor at different nucleotides during Drosha and Dicer processing. Another reason of this variability could be the nucleotide replacement and additions or deletions in the 59 and 39 miRNA ends. All this posttranscriptional mechanisms, called RNA editing, generate the isomiRs, defined as sequence variations in mature miRNAs and commonly reported in miRNA studies [10,70]. Mechanisms of miRNA sequence diversification are extensively reviewed by Berezikov [16], suggesting a functional meaning to all these processes. Furthermore, different patterns described on isomiR expression suggest different strategies in the regulation of the expression. Assuming the expression data obtained in this study, there are miRNAs where the main isomiR could represent the major regulation function and there are miRNAs which their regulation function could be assumed by more than one isomiR (Table S3).
Considering the proportion of unique sequences obtained per breed in reference to the total sequences, Iberian and Vietnamese breeds obtained the highest proportion (38% and 39%, re- spectively) showing that, if the number of miRNAs was not increased in these libraries, what was increased was the number of described isomiRs. The percentage of unique sequences in the other libraries and also in the whole data (all libraries) was between 15 and 19%. Thus, the variability in miRNA sequences seemed to be higher in European breeds like Iberian breed and in Asian breeds such as Vietnamese breed. In contrast, when European commercial breeds were created through the European and Asian breeds introgression, their variability could be reduced due to their strong selection in production traits, as it is indicated if we consider the number of isomiRs in each European commercial breed studied.
Looking at the most expressed isomiR by miRNA in each breed we can determine that largely it was the same isomiR in all studied breeds. However, in 76 miRNAs, like Ssc-miR-99a, the most expressed isomiR varies depending on the breed, being more expressed one isomiR (U0009355: 59-AACCCGTAGATCC-GATCTTGTG-39) in European and Asian breeds, and other isomiR (U0009353: 59-AACCCGTAGATCCGATCTTGTGA-39) in European commercial breeds. In other cases, there is an agreement in all breeds for the most expressed isomiR except in one breed, like in Ssc-miR-21 and Ssc-miR-100 in Iberian breed, Ssc-miR-152 in Piétrain breed, Has-miR-200c-3p in Vietnamese breed or Hsa-let-7d-5p in Landrace breed. There are also few cases where the most expressed isomiR is different in almost each breed, such as in Ssc-miR-199a* and Ssc-miR-423-5p. All these different situations suggest that the genetic variability degree between breeds is high and it could play an important role in the post-transcriptional regulation mechanisms leading to perform different functions.
Another variation factor was that there were some differences between the most expressed isomiR in each miRNA and its concordance with the described sequence in miRBase database, diverging them in 72 miRNAs found in this study. As example, the most representative sequence for Ssc-miR-23a in this study was 59-ATCACATTGCCAGGGATTTCCA-39, found 3,116 times and it corresponded to Bta-miR-23a. The sequence described in miRBase database for Ssc-miR-23a (59-ATCACATTGCCAGG-GATTTCC-39) was found only 622 times (Table 4). It also suggests that there is not a fixed predominant isomiR, but there is variability in isomiR expression given not only by many factors like age, tissue or disease [71], but also by other factors like breed [50].
Despite of all miRNA homologies found, this study was also focused on describing novel miRNAs. The strict methodology used proposed 8 novel miRNAs to be validated through RT-qPCR and finally 5 miRNAs were successfully validated. These confirmed novel porcine miRNAs were mapped once in the pig genome (Table S7) and their pre-miRNA folding were successfully predicted ( Figure 2). Moreover, three out of the thirteen selected miRNAs from the kidney miRNAome to be amplified through RT-qPCR were orthologous (Hsa-miR-200b-3p, Hsa-miR-200c-3p and Bta-miR-193b, Table 5) and, therefore, they were also confirmed as new porcine miRNAs. Significant differential expression regarding breed groups was obtained by RT-qPCR in eight miRNAs: Hsa-miR-200b-3p, Hsa-miR-200c-3p, Ssc-miR-126, Ssc-miR-126*, Ssc-miR-99a, Bta-miR-193b, Ssc-miR-486 and Ssc-let-7f. However, expression results from RT-qPCR were not in agreement with the differential expression study from high throughput sequencing data in some cases. These differences could be explained by many factors, like some possible bias, the fact that not each miRNA is amplified in the same proportion with high throughput sequencing [72,73] and the presence of isomiRs, which make difficult to perform an accurate filtering of sequences. Other factors to be considered are some diversity and complexity in the miRNAs nomenclature between species listed in the miRNAs databases [15], which make difficult to perform a proper sequence classification. In this sense, cluster analysis methodology for novel miRNA detection was chosen to avoid problems with miRNA nomenclature. It is also important to consider the specificity in qPCR amplification, where primers designed in each miRNA were exclusive for the most expressed isomiR in each miRNA, amplifying only that one. Overall, the relevant topic of next generation sequencing about if it can be treated as expression data taking into account its large amount of reads, is still discussed [74]. The slightly different expression results between high throughput sequencing and RT-qPCR data reveal a change of expression tendency in some miRNAs evidencing the ''semiquantitative'' nature of high throughput sequencing methods. Low correlations (from 0.52 to 0.003) were determined between high throughput sequencing by 454 (Roche, Germany) and RT-qPCR data. In this sense, the bias in this study could be explained by the lower number of sequences obtained. For instance, Hsa-miR-200b-3p appears in a higher expression in Asian breeds in sequencing data while it is more expressed in European breeds according to real time RT-qPCR data. Furthermore, Ssc-miR-126* was more represented in European commercial breeds, whereas in RT-qPCR it was up regulated in European breeds. The same case happened with Ssc-miR-126, a conclusive fact taking into account that both miRNAs are transcribed together.
There are similar cases with low expressed miRNAs (according to sequence count) such as Ssc-let-7f, which was more expressed in Asian breeds but appears up regulated in European breeds in real time RT-qPCR data.
Analysing the target pathways of differentially expressed miRNAs between breed groups, it was found that Hsa-mir-200b-3p, the most expressed miRNA in this study and up regulated in European breeds, was related in several kidney diseases like tubulointerstitial fibrosis [75] or hypertensive glomerulosclerosis [12], elucidating the importance of this miRNA in kidney physiology. Targets of Hsa-miR-200b-3p were involved in several pathways like renal cell carcinoma, associated to some oncogenes such as MET, or tumors suppressors like VHL, FH and BHD. Interestingly, miRNA targets were related with the vascular smooth muscle contraction and calcium signaling pathways. Both pathways are important in muscular growth and development in which Hsa-miR-200b-3p, being more expressed in European pig breeds, could contribute to a greater efficiency. Furthermore, European commercial breeds are precisely selected for production traits and strongly characterized for animal production. Kidney and muscle are strongly related because kidney is in charge of secreting erythropoietin, the hormone responsible to activate erythropoiesis and also to transport oxygen to the muscles and For HTS data, fold changes from sequence counting between breed groups were calculated from normalised data in counts per thousand for each library and averaged per groups. Positive and negative signs indicate that the level of gene expression is higher for the first or the second group of the test, respectively. Analysis of variance (ANOVA) including breed as fixed factor was performed. Significance was set at P,0.05. Fold changes for clusters could not be calculated due to their low total counts. For qPCR data expression study, the quantity obtained of each miRNA in each sample was normalised by the Normalization Factor and corrected in relation to the lowest normalised value. Analysis of variance (ANOVA) including the RT and breed as fixed factors was performed and fold change from Least Squares Means (LSM) between breed groups were calculated. Significance was set at P,0.05. Scheffe test determined whether there was significant differential expression between breed groups ( 1 : suggestive p-value,0.1, * : significant with p-value,0.05, ** : significant with p-value,0.01, *** : significant with p-value,0.001). doi:10.1371/journal.pone.0055402.t006 other tissues, and thus, favoring muscular growth. Hsa-miR-200c-3p, which belongs to the same family that Hsa-miR-200b-3p, is also differentially expressed being up regulated in Asian breeds. Some studies have related it to some mechanisms of cancer development in many organs, such as pancreas, bladder, ovaries, prostate or breast. [76][77][78][79][80]. Functional analysis related targets of Hsa-miR-200c-3p with neoplasias, specifically in renal cell carcinoma. Furthermore, Hsa-miR-200c-3p is also involved in some pathways related to reproduction, like oocyte meiosis or progesterone-mediated oocyte maturation. Breeds with an Asian origin have better phenotypes for reproduction traits and, interestingly, Asian breeds have Hsa-miR-200c-3p up regulated suggesting its involvement in the regulation of the reproduction traits. However, Ssc-miR-99a and Hsa-miR-200b-3p targets are also related to these reproduction pathways, both up regulated miRNAs in European breeds which could play an antagonistic role in European breeds. Some targets of Ssc-miR-126, a down regulated miRNA in Asian breeds, participate in cellular processes like focal adhesion and regulation of actin cytoskeleton, while Ssc-miR-126*, also down regulated in Asian breeds, regulate targets implicated with cellular processes like adherens, gap and tight junction. Moreover, both miRNAs also have a strong relation with many cancer pathways, which were related in all miRNAs. Ssc-let-7f is an ubiquitous miRNA from let-7 family, related to many pathways involved in wide range of physiological functions (cellular, metabolic and environmental information processes) and also to many diseases, like renal cell carcinoma [81,82].
This study opens a wide field about the miRNA specificity, which is not only by tissue, age or species level, but also breed factor is crucial to manage the phenotypic changes and, therefore, the genetic variability. The role of miRNAs in gene expression regulation is much more complex than it was thought initially, being necessary to study in depth which roles isomiRs play and to uncover all factors that make the post-transcriptional regulation through miRNAs so complex.

Supporting Information
Table S1 Adaptors used for the construction of each library. 1 : Forward adaptor, 2 : Reverse adaptor. The 5 nt code used for each library is in bold in the adaptor sequence. Reverse adaptor was used in all breeds libraries. (DOC)