Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 https://doi.org/10.1186/s12711-021-00632-3 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access Variability in porcine microRNA genes and its association with mRNA expression and lipid phenotypes Emilio Mármol‑Sánchez1, María Gracia Luigi‑Sierra1, Anna Castelló1,4, Dailu Guan1, Raquel Quintanilla2, Raul Tonda3 and Marcel Amills1,4* Abstract Background: Mature microRNAs (miRNAs) play an important role in repressing the expression of a wide range of mRNAs. The presence of polymorphic sites in miRNA genes and their corresponding 3′UTR binding sites can disrupt canonical conserved miRNA–mRNA pairings, and thus modify gene expression patterns. However, to date such poly‑ morphic sites in miRNA genes and their association with gene expression phenotypes and complex traits are poorly characterized in pigs. Results: By analyzing whole-genome sequences from 120 pigs and wild boars from Europe and Asia, we identified 285 single nucleotide polymorphisms (SNPs) that map to miRNA loci, and 109,724 SNPs that are located in predicted 7mer-m8 miRNA binding sites within porcine 3′UTR. In porcine miRNA genes, SNP density is reduced compared with their flanking non-miRNA regions. By sequencing the genomes of five Duroc boars, we identified 12 miRNA SNPs that were subsequently genotyped in their offspring (N = 345, Lipgen population). Association analyses of miRNA SNPs with 38 lipid-related traits and hepatic and muscle microarray expression phenotypes recorded in the Lipgen population were performed. The most relevant detected association was between the genotype of the rs319154814 (G/A) SNP located in the apical loop of the ssc-miR-326 hairpin precursor and PPP1CC mRNA levels in the liver (q-value = 0.058). This result was subsequently confirmed by qPCR (P-value = 0.027). The rs319154814 (G/A) genotype was also associated with several fatty acid composition traits. Conclusions: Our findings show a reduced variability of porcine miRNA genes, which is consistent with strong purifying selection, particularly in the seed region that plays a critical role in miRNA binding. Although it is gener‑ ally assumed that SNPs mapping to the seed region are those with the most pronounced consequences on mRNA expression, we show that a SNP mapping to the apical region of ssc-miR-326 is significantly associated with hepatic mRNA levels of the PPP1CC gene, one of its predicted targets. Although experimental confirmation of such an interac‑ tion is reported in humans but not in pigs, this result highlights the need to further investigate the functional effects of miRNA polymorphisms that are located outside the seed region on gene expression in pigs. Background Mature microRNA transcripts (miRNAs) are short (~ 22 nt) non-coding RNAs that play an essential role in the regulation of gene expression [1]. During the biogen- *Correspondence: marcel.amills@uab.cat 1 esis of miRNAs, one strand of the miRNA duplex binds Centre for Research in Agricultural Genomics (CRAG), CSIC‑IRTA‑UAB‑UB, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain to the guide-strand channel of an Argonaute protein Full list of author information is available at the end of the article forming a miRNA-induced silencing complex (miRISC), © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​ mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​mmons.​org/​publi​cdoma​in/​ zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 2 of 21 which has the ability to repress mRNA expression by their maturation and loading into the functional silencing binding to specific 3′UTR target sites [2]. In postem- complex [18, 19]. Moreover, Saunders et al. [5] showed bryonic cells, this repressor mechanism acts mainly by that a broad array of predicted miRNA target sites in the destabilizing the mRNA through decapping and poly(A)- 3′UTR of mRNAs are polymorphic, which suggests that tail shortening [3, 4] and less often by hindering transla- purifying selection in these regions is less intense than in tion [2]. The binding of the miRNA to its 3′UTR target miRNA genes [17]. Overall, these findings provide evi- site depends critically on the sequence of the seed region, dence for the existence of functional conserved mecha- which encompasses nucleotides (nt) 2–8 from the 5′ nisms of miRNA-mediated gene regulation that are end of the mature miRNA and interacts with the target influenced by polymorphic sites at both miRNA loci and site through Watson–Crick pairing [1]. Polymorphisms their 3′UTR binding sites [20]. located within the seed region can affect the proper rec- The wild ancestors of pigs were independently domes- ognition of mRNA targets [5, 6]. Nevertheless, imper- ticated in the Near East and China ~ 10,000 years before fect seed matches can be compensated by nucleotides 13 present [21]. Then, domestic pigs spread worldwide to 16 of the mature miRNA, thus providing additional to become one of the most important sources of ani- anchor pairing to the seed region [2]. Other sites that are mal protein for humans and they have diversified into relevant for miRNA processing and function are the basal an extensive array of breeds with distinct morphologi- UG, flanking CNNC and apical UGU motifs [7], as well cal and productive features [22]. Pig phenotypes might as a mismatched GHG motif [8], all of which contrib- be explained, at least partially, by modifications in the ute to facilitate miRNA processing [2]. Variability in the microRNA-mediated regulation of gene expression 3′UTR miRNA binding sites can also affect gene expres- [23]. Indeed, several studies have reported associations sion at the post-transcriptional level [5, 9, 10]. between SNPs in miRNA binding sites and porcine Evolutionary conservation of miRNAs across species phenotypic variation [24–29], while fewer studies have is associated with their level of expression [11, 12] and investigated the association between SNPs in miRNA with the functional importance of the regulatory net- genes and complex traits in pigs [30–33]. For instance, works that they modulate [13, 14]. In addition, conser- Chai et al. [33] reported a significant causal relationship vation of the repertoire of mRNAs that are targeted by between a SNP in the seed of ssc-miR-378-3p and exten- a given miRNA depends mostly on the age of the asso- sive modifications of its target mRNA repertoire. Moreo- ciated miRNA gene, with novel miRNAs acquiring tar- ver, they associated the presence of the observed miRNA gets more rapidly than ancient miRNAs [15]. Conserved seed variant to a structural reconfiguration of the hairpin, sites are particularly enriched in the 5′end of 3′UTR [12], leading to an enhanced expression of the mature miRNA although selective pressure on miRNA binding sites is [33]. not uniform across miRNA genes [16]. The polymorphism of miRNA genes has not been sys- Saunders et al. [5] investigated the variability of 474 tematically characterized in pigs despite its potential human miRNA genes and found that the density of sin- impact on gene regulation and phenotypic variation. In gle nucleotide polymorphisms (SNPs) is lower in miRNA the current work, our aim was to elucidate the patterns loci than in their flanking regions. Moreover, Saunders of the variability of miRNA genes in European and Asian et al. [5] found that ~ 90% of the human miRNA genes wild boars and domestic pigs and to infer whether such do not contain polymorphisms, and that most of the patterns are influenced by purifying selection. We also SNPs mapping to miRNAs are located outside the seed investigated whether SNPs in miRNA genes are associ- region, which provides evidence that the variability of ated with liver and muscle gene expression traits and this critical functional motif evolves under strong selec- lipid phenotypes recorded in a Duroc pig population. tive constraints [5]. Indeed, polymorphisms in the first 14 nucleotides of mature miRNAs, and particularly those Methods within the seed region, might abrogate the binding of Characterization of polymorphisms in miRNA genes miRNAs to their 3′UTR targets, leading to an extensive and their 3′UTR binding sites in a worldwide sample rewiring of miRNA-mediated regulatory networks [17]. of pigs and wild boars Mammalian miRNA knockouts often display abnormal Retrieval of porcine whole‑genome sequences phenotypes, reduced viability, and clinical disorders [2], Whole-genome sequences (WGS) from 120 wild boars although functional redundancy among miRNA fam- and domestic pigs (Sus scrofa) were retrieved from the ily members might mitigate, to some extent, the sever- NCBI Sequence Read Archive (SRA) database (https://​ ity of such manifestations [2]. Polymorphisms within www.​ncbi.​nlm.​nih.​gov/​sra). Most of these WGS have miRNA loci that lie outside the seed region can also been reported in previous publications [34–37] and affect the processing and stability of miRNAs during detailed information is available in Additional file 1: Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 3 of 21 Table S1. The 120 WGS corresponded to Asian domes- reliability and reproducibility of the detection of InDels tic pigs (ADM, N = 40), Asian wild boars (AWB, N = 20), based on commonly used pipelines for variant calling are European domestic pigs (EDM, N = 40) and European low across tools and analyzed datasets [42, 43]. wild boars (EWB, N = 20). The ADM group included WGS from Meishan, Tongchen, Jinhua, Rongchan, Identification of single nucleotide polymorphisms in miRNA Wuzhishan, Tibetan, Sichuan, Hetao, Minzhu, Bamaix- genes and in their 3′UTR binding sites ang and Laiwu pigs, and the EDM group was composed Single nucleotide polymorphisms mapping to annotated of WGS from Pietrain, Mangalitza, Iberian, Duroc, porcine miRNA loci (N = 370) were retrieved by using American Yucatan (from America, but with a European the curated Sscrofa11.1 annotation for miRNA regions origin), Yorkshire, Landrace, Hampshire and Large-white available in the miRCarta v1.1 [44] and miRBase [45] pigs. Both AWB and EWB samples were selected accord- databases. In addition, annotated mature miRNA loci ing to their geographic distribution (see Additional file 1: (N = 409) within miRNA genes (N = 370) were retrieved, Table S1). The EWB group included one sample from the and SNPs located in the mature and seed regions (posi- Near East. Raw data in SRA format were downloaded tions 2 to 8 at the 5′ end of the mature miRNA) were from SRA public repositories and converted into fastq identified. A comprehensive description of miRNA gene format by using the fastq-dump 2.8.2 tool that is avail- regions and miRNA–mRNA interactions is shown in able in the SRA-toolkit package (https://​ncbi.​github.​io/​ Fig. 1a, b. sra-​tools/). We also retrieved the annotated 3′UTRs from por- cine mRNA transcripts available at Ensembl reposi- tories (v.100, http://​www.​ensem​bl.​org/​info/​data/​ftp/​ Whole‑genome sequence data processing and calling index.​html) and the corresponding set of sequences of single nucleotide polymorphisms was interrogated against the seed regions of the 409 Fastq paired-end files generated from SRA data were porcine mature miRNAs considered in our study. Seed filtered according to their quality, and sequence adapt- sequences were reverse-complemented and searched ers were trimmed by using the Trimmomatic software along the 3′UTR sequences of mRNA genes by using (v.0.36) with default parameters [38]. Trimmed paired- the locate function of the SeqKit toolkit [46]. All 7mer- end sequences were aligned against the Sus scrofa refer- m8 canonical seed pairing sites and the 8mer seed pair- ence genome (Sscrofa11.1) [39] with the BWA-MEM ing sites (Fig. 1) were identified for each miRNA seed, algorithm [40] and default settings. Sequence alignment and, based on this information, the corresponding puta- map (SAM) formatted files were sorted and transformed tive miRNA–mRNA target pairs were established. Addi- into binary (BAM) formatted files. PCR duplicates were tional non-canonical miRNA–mRNA interactions were subsequently identified with the MarkDuplicates package not considered (Fig. 1b) as they are regarded as less bio- from the Picard tool (https://​broad​insti​tute.​github.​io/​pic- logically significant than 7mer-m8 or 8mer interactions ard/) and removed to perform insertion-deletion (InDel) [2, 47, 48]. Subsequently, the genomic locations of both realignment with the IndelRealigner package from the 7mer-m8 and 8mer matching regions were determined, Genome Analysis Toolkit (GATK v.3.8) [41]. Base qual- and SNPs that were located in such predicted miRNA ity score recalibration (BQSR) was implemented with the binding sites were identified. BaseRecalibrator package from GATK v.3.8. Variant call- Principal components analyses (PCA) were performed ing was implemented with the HaplotypeCaller function with the smartPCA software [49]. These PCA analy- according to GATK best practices [41]. Individual gVCF ses were based on: (1) autosomal whole-genome SNPs, formatted files, including both polymorphic and homozy- (2) autosomal SNPs from miRNA genes, (3) autosomal gous blocks, were generated with the GenotypeGVCFs SNPs from 3′UTRs, or (4) autosomal SNPs from 8mer package from the GATK v.3.8 tool [41], and they were and 7mer-m8 3′UTR sites. When implementing the PCA merged into separate multi-individual variant calling for- based on whole-genome data, SNPs that complied with mat (VCF) files containing single polymorphic and InDel the following parameters were retained: (1) minimum sites for each defined population (i.e. ADM, EDM, AWB allele frequency > 0.05, and (2) Hardy–Weinberg equi- and EWB). Finally, variants were distributed into SNP librium exact test P-value > 0.001. In contrast, no filter- and InDel files with the SelectVariants package and qual- ing was performed in the PCA based on miRNA SNPs ity-filtered using the FilterVariants function from GATK because of their small number. v.3.8 [41]. For SNPs, we used the following parameters: QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < − 12.5, Frequency and distribution of SNPs in miRNA genes ReadPosRankSum < − 8.0, SOR > 3.0. In the current study, To investigate the distribution of SNPs along miRNA we did not take InDels into consideration because the genes, SNPs were classified according to their location, Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 4 of 21 Fig. 1 Schematic representations of a precursor miRNA sequence hairpin conformation, with mature miRNAs located in the stems (higlighted in blue) and their corresponding seed regions (in red). b Canonical miRNA–mRNA interactions where each binding site is conformed by six to seven Watson–Crick pairings (vertical lines) to the seed region of the mature miRNA (2nd to 8th 5′ nucleotides). As reported elsewhere [2, 47, 48], 7mer-m8 and 8mer binding sites are involved in the most functionally active canonical miRNA–mRNA interactions i.e., SNPs within the seed region of mature miRNAs were to obtain the mature length (6053 bp), which corresponds flagged as “seed”, whereas SNPs located in mature miR- to the total length of all mature miRNA sequences exclud- NAs but outside the seed were classified as “mature”. In ing the seeds. The precursor length (22,229 bp) was calcu- addition, SNPs classified as “mature” were assigned to lated by subtracting the seed length (2863 bp) and mature the following subtypes: (1) “anchor” (1st position at the length (6053 bp) from the total length (31,145 bp) of all 5′ end), and (2) “supplemental pairing” (position 13 to 18 miRNA loci (N = 370). Then, the SNP density ( D) in each from the 5′ end). The remaining SNPs were assigned to of these regions was calculated as follows: the “precursor” class. Allele frequencies were estimated in Nr × 100 the whole set of porcine samples (N = 120) and in each D= , EDM, ADM, EWB and AWB groups independently (see Lr Additional file 2: Table S2). where Nr is the total number of miRNA SNPs ( N = 285), In order to calculate the SNP density in precursor, or the number of SNPs detected in each of the defined mature and seed miRNA regions, first we calculated the regions, (i.e., 221, 52 and 12 SNPs in precursor, mature total nucleotide length of each of these regions in the set and seed regions, respectively); Lr is the total nucleo- of analyzed porcine miRNA loci. Seed length was cal- tide length of all miRNA loci (31,145 bp), or the total culated by considering seven nucleotide positions in number of nucleotides in miRNA loci belonging to the each of the annotated mature miRNAs (N = 409) map- “precursor” (22,229 bp), “mature” (6053 bp) or “seed” ping to the 370 miRNA loci considered in our study (2863 bp) regions. These calculations yield the num- (7 bp × 409 = 2863 bp). The total length of the 409 mature ber of SNPs per bp for each category, but we decided to miRNAs was determined (~ 22 bp × 409 = 8916 bp), and adjust this estimate to a 100-bp window, which results the total seed length (2863 bp) was subtracted from 8916 bp Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 5 of 21 in the multiplication of Nr by 100 in the above formula. piglets were weaned, castrated and subsequently fattened We also compared the SNP density between miRNA loci in four contemporary batches on the IRTA pig experi- ( N = 370) and their flanking genomic regions (± 1 kb) mental farm in Monells (Girona, Spain) under intensive by considering 10 upstream and downstream nucleotide standard conditions. Once they reached ~ 122 kg of live bins of 100 bp each. Differences in SNP density across weight (~ 190 days of age), they were slaughtered in a miRNA loci, as well as its comparison with surrounding commercial abattoir following the guidelines of the cur- regions, were estimated with a non-parametric approach rent Spanish legislation (https://​www.​boe.​es/​buscar/​ using a Mann Whitney U test [50]. doc.​php?​id=​BOE-A-​1995-​3942). After slaughtering, tis- sue samples from the gluteus medius (GM) and longissi- Investigating the association of miRNA polymorphisms mus dorsi (LD) skeletal muscles and liver were obtained with gene expression and phenotype data recorded according to previously described methods [54–56]. in Duroc pigs Total DNA was extracted from each sample following Whole‑genome sequencing of five Duroc pigs Vidal et al. [51]. In total, 345 DNA samples from the ini- In 2003, five Duroc boars were selected and used as tial set of 350 pigs were successfully obtained and pro- founders of a half-sib population of purebred Duroc cessed for further genotyping. Total RNA was extracted pigs that are dedicated to the production of high-quality from GM (N = 89 pigs) and liver (N = 87 pigs) tissue sam- cured ham. By whole-genome sequencing of these five ples by using the Ribopure RNA isolation kit (Ambion, individuals, our aim was to identify segregating miRNA Austin, TX). Expression mRNA profiles were charac- SNPs and investigate their association with gene expres- terized by hybridizing total RNA to GeneChip porcine sion traits and lipid phenotypes recorded in the offspring arrays (Affymetrix Inc., Santa Clara, CA), which encom- of the five boars (Lipgen population, N = 350). Genomic pass 23,998 probes [54, 57]. Further details about tissue DNA was extracted [51] and sequenced at the Centre collection, sample selection, RNA isolation and microar- Nacional d’Anàlisi Genòmica (CNAG, Barcelona, Spain). ray hybridization procedures are in Cánovas et al. [54]. Paired-end multiplex libraries were prepared according Microarray data pre-processing, background correc- to the instructions of the manufacturer with the KAPA tion, normalization and ­log2-transformation of expres- PE Library Preparation kit (Kapa Biosystems, Wilming- sion estimates were performed with a robust multi-array ton, MA). Libraries were loaded to Illumina flow-cells average method per probe [58]. The mas5calls function for cluster generation prior to producing 100-bp paired- from the affy R package [58] was then applied to iden- end reads on a HiSeq2000 instrument following the Illu- tify probes with intensity signals that were higher than mina protocol. Base calling and quality control analyses the background noise. This function applies a Wilcoxon were performed with the Illumina RTA sequence analysis signed rank-based gene expression presence/absence pipeline according to the instructions of the manufac- detection algorithm for labeling expressed probes in turer. Quality-checked filtered reads were mapped to the each sample. Control probes and those with expression Sus scrofa genome version 11.1 and processed for SNP levels below the detection threshold in more than 50% calling according to GATK best practice recommenda- of samples were discarded from further analyses. As a tions [41] and following the same protocol described in result, we retained only the probes which corresponded the previous section entitled “Whole-genome sequence to mRNAs that were robustly expressed in GM muscle data processing and calling of single nucleotide poly- (12,131 probes) and liver (12,787 probes) tissues. Probes morphisms”. Since only five boars were sequenced, SNPs from the GeneChip Porcine Genomic array identifiers within miRNA loci were retrieved without applying qual- (Affymetrix, Inc., Santa Clara, CA) of GM skeletal mus- ity filtering criteria, in order to maximize the number of cle and liver tissues were assigned to mRNA genes anno- identified miRNA SNPs. tated in the Sscrofa11.1 assembly [39] using the BioMart tool [59]. Recording of phenotypes related with gene expression With regard to lipid-related phenotypes, 38 traits were and lipid traits in the Lipgen population measured in the Duroc Lipgen population, i.e. backfat In total, 350 Duroc barrows, sired by the five Duroc thickness measured between the 3rd and 4th ribs, backfat founder boars mentioned above, were used as a resource thickness at the last rib, ham fat thickness, and intramus- population (Lipgen population [52, 53]) to evaluate asso- cular fat content and composition of GM and LD skeletal ciations between miRNA SNPs and phenotypes related muscle samples (N = 345) [53, 55] (see Additional file 3: with gene expression and lipid traits. The five sequenced Table S3). Briefly, intramuscular fatty acid (IMF) content boars were mated with 400 sows on three farms and one in the GM and LD muscles was measured with the Near offspring per litter was selected for phenotypic recording Infrared Transmittance technique (NIT, Infratec 1625, (only 350 individuals provided valid records). All selected Tecator Hoganas, Sweden), whereas muscle cholesterol Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 6 of 21 measurements were inferred following Cayuela et al. [60]. PLINK software [61] in order to generate formatted files Gas chromatography of methyl esters was used to deter- for subsequent analyses. The Genome-Wide Efficient mine muscle fatty acids (FA) composition, including the Mixed-Model Association (GEMMA) software [62] was percentages of saturated (SFA), unsaturated (UFA), mon- used to implement association analyses between geno- ounsaturated (MUFA) and polyunsaturated (PUFA) fatty typed SNPs and phenotypes related with microarray gene acids. Live and carcass weights (used as covariates in the expression data and lipid traits. The following univariate statistical analyses), as well as backfat and ham fat thick- mixed model was used: ness, were measured on a regular basis, prior (live weight y = Wα + xδ + u + ε, and backfat thickness) and after slaughter (carcass weight and ham fat thickness). Mean and standard deviations of where y is the vector of recorded phenotypes for each lipid traits recorded in the Lipgen population (N = 345) individual; α is a vector indicating the intercept plus the are described in Additional file 3: Table S3. considered fixed effects, i.e., batch effect with four cat- egories (all traits), farm of origin effect with three cat- Genotyping of a panel of single nucleotide polymorphisms egories (all traits) and laboratory of processing with two mapping to microRNA genes in the Lipgen population categories (GM and liver microarray expression data). α Twelve SNPs that mapped to miRNA loci and segre- also includes the following covariates: IMF of LD (for LD gated in the five sequenced parental Duroc boars were fatty acid composition traits), IMF of GM (for GM fatty genotyped in the 345 Duroc pigs from the Lipgen popu- acid composition traits), live weight (for backfat thick- lation (Table 1). Briefly, selected miRNA SNPs and their ness) and carcass weight (for ham fat thickness); W is the flanking regions (60 bp upstream and downstream) were incidence matrix relating phenotypes with their corre- used for the assay design with the Custom TaqMan Assay sponding effects; x is the vector of genotypes for the 12 Design Tool website (https://​www5.​appli​edbio​syste​ms.​ miRNA polymorphisms; δ is the allele substitution effect com/​tools/​cadt/; Life Technologies). Genotyping was for each polymorphism; u is a vector of random indi- carried out at the Servei Veterinari de Genètica Molecu- vidual genetic effects with an n-dimensional  multivariate lar of the Universitat Autònoma of Barcelona (http://​sct.​ normal distribution MVNn ∼ 0, τ −1 K , where τ −1 is uab.​cat/​svgm/​en) by using a QuantStudio 12K Flex Real- the variance of the residual errors,  is the ratio between Time PCR System (Thermo Fisher Scientific, Barcelona, the two variance components and K is the known relat- Spain). edness matrix derived from SNP information; and ε is the vector of residual errors. Association analyses between miRNA SNPs and mRNA Association analyses were performed between each expression and lipid phenotypes miRNA polymorphism (Table 1) and the 38 lipid-related Genotype data corresponding to the 12 miRNA SNPs traits listed in Additional file 3: Table S3. Moreover, we mentioned above (Table 1) were processed with the explored the associations between miRNA SNPs and Table 1 List of miRNA polymorphisms genotyped in a population of Duroc pigs (N = 345, Lipgen population) microRNA SSCa Start (bp) End (bp) Strand SNP Type Alt. ­alleleb Frequencyc ssc-miR-339 1 164,025,972 164,026,057 − rs81349391 Apical loop G 0.4138 ssc-miR-130a 2 13,296,695 13,296,773 − rs344472188 Apical loop G 0.1364 ssc-miR-23a 2 65,308,117 65,308,186 + rs333787816 Precursor stem C 0.5398 ssc-miR-30d 4 6,948,669 6,948,747 + rs340704946 Precursor stem G 0.4762 ssc-miR-371 6 56,427,208 56,427,285 − rs320008166 Precursor stem C 0.2955 ssc-miR-429 6 63,491,921 63,492,001 + rs323906663 Precursor stem A 0.2955 ssc-miR-326 9 9,581,944 9,582,034 − rs319154814 Apical loop A 0.5852 ssc-miR-34c 9 39,280,278 39,280,357 + rs321151601 Mature region A 0.0805 ssc-miR-378–2 12 36,947,443 36,947,510 + rs341950320 Precursor stem A 0.1176 ssc-miR-15b 13 100,083,172 100,083,269 + rs334680106 Precursor stem T 0.2706 ssc-miR-486 17 10,758,818 10,758,899 − rs335924546 Precursor stem T 0.3391 ssc-miR-335 18 18,341,568 18,341,659 − rs334590580 Precursor stem C 0.1875 a SSC: porcine chromosome b Alt. allele: alternative allele of the genotyped SNP c Frequency: frequency of the alternative allele in the Lipgen population Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 7 of 21 expression levels of mRNAs that fulfilled the follow- software [67] was used for determining enriched Reac- ing conditions: (1) the 3′UTR of the mRNA contains at tome and KEGG pathways. A one-sided hypergeometric least one 7mer-m8 site complementary to the seed of test of significance was applied for determining enriched the miRNA harboring the SNP to be tested (as predicted terms and multiple testing correction was implemented with the locate tool from SeqKit software [46]), and (2) with the FDR method [65]. the existence of experimentally-validated mRNA-miRNA interactions between the given mRNA and the tested Confirmation of associations between ssc‑miR‑326 miRNA have been confirmed in humans according to rs319154814 genotypes and gene expression data information provided in the DIANA-Tarbase v8 database by RT‑qPCR [63]. The first condition was established because we were The hepatic mRNA levels of three of the mRNA tran- interested in focusing our analyses on miRNA–mRNA scripts among the most significantly associated with pairs for which a direct physical interaction (through ssc-miR-326 rs319154814 (G/A) genotype were analyzed complementarity between the seed of the miRNA and by reverse transcription-quantitative polymerase chain a target site in the 3′UTR of the mRNA) is feasible. The reaction (RT-qPCR). The β-actin (ACTB) and TATA- second condition was imposed because there is exten- box binding protein (TBP) genes were used as endog- sive evidence that in silico prediction or miRNA binding enous controls, as previously reported [68]. In brief, sites can have high false positive rates [64]. We consid- 100 mg of liver tissue from 10 selected samples (5 from ered as valid miRNA–mRNA interactions those that have each GG and AA genotypes) were placed in 1 mL TRI- been experimentally validated in humans by using wet zol reagent (Invitrogen Corp., Carlsbad, CA, USA) and lab methods such as cross-linking, ligation and sequenc- subsequently homogenized using the Lysing Matrix D ing of hybrids (CLASH), photoactivatable ribonucleo- reagent (MP Biomedicals, Santa Ana, CA) in a Precellys side-enhanced and high-throughput sequencing of RNA 24 tissue homogenizer (Bertin Instruments, Rockville, isolated by crosslinking and immunoprecipitation (PAR- MD). Total RNA isolation was performed following the CLIP and HITS-CLIP) and luciferase assays [63]. protocol described by Rio et al. [69]. Reverse transcrip- It should be noted that the seeds of the porcine miR- tion was achieved with the High-Capacity cDNA Reverse NAs harboring the 12 genotyped SNPs were completely Transcription Kit (Applied Biosystems, Foster City, conserved with regard to the seeds of the human ortholo- CA) following the instructions of the manufacturer. We gous miRNA genes (see Additional file 4: Figure S1). This used 1 μg of total RNA as a template in a final volume of feature further supports the extrapolation of miRNA– 20 μL and the synthesized cDNA was diluted to 1:20 in mRNA interactions experimentally validated in humans Milli-Q water. One pair of primers spanning exon-exon [63] to pigs. For the sake of completeness, we also car- junctions was designed for each gene (see Additional ried out a complementary association analysis between file 5: Table S4) with the Primer Express software (Life miRNA SNPs and the mRNA levels of the whole sets of Technologies Corporation). Real-time qPCR reactions genes expressed in the GM and liver tissues (i.e., without contained 10 μL of SYBR Select Master Mix, 300 nM of applying the two conditions mentioned above) but the each primer, and 5 μL of 1:20 cDNA dilution, in a final results of such analyses should be interpreted with cau- volume of 20 μL. Reactions were run in an ABI PRISM tion due to reasons stated above. 7900HT instrument (Applied Biosystems, Foster City, In association analyses with lipid traits or gene expres- CA). The thermal cycle was set as follows: 2 min at 50 °C, sion phenotypes, the alternative hypothesis H ­ 1: δ ≠ 0 one denaturing step at 95 °C during 10 min, followed by was contrasted against the null hypothesis H ­ 0: δ = 0 with 40 cycles of 15 s at 95 °C and 1 min at 60 °C. Moreover, a a likelihood ratio test. The statistical significance of the melting curve analysis was performed (i.e., 95 °C for 15 s, associations between miRNA SNPs and lipid and mRNA 60 °C for 15 s and a gradual increase in temperature, with expression phenotypes was assessed with a false discov- a ramp rate of 1% up to 95 °C, followed by a final step ery rate (FDR) approach [65]. of 95 °C for 15 s), in order to assess the specificity of the reactions. A standard curve with serial dilutions from a pool including all the cDNA samples was used to assess Pathway enrichment analysis whether amplification efficiencies of the three genes The lists of probes/genes that are significantly associ- under analysis were comprised in the 90 to 110% range. ated at the nominal level (P-value < 0.05) with miRNA All reactions were run in triplicate. SNPs genotypes (N = 12) after considering the criteria for The ­2−ΔΔCt method [70] was used for the relative miRNA target inference (see “Methods” above) were used quantification of gene expression (Rq), using the group as inputs for pathway enrichment analyses. The ClueGO of GG samples as a calibrator. Subsequently, Rq values v2.5.0 plug-in application [66] within the Cytoscape 3.8.2 were ­log2 transformed and the significance of expression Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 8 of 21 differences between GG and AA genotypes was assessed is available in Additional file 6: Table S5. Most of these using the Welch’s t-test for unpaired groups of samples SNPs were biallelic (96.82%), 777,008 (1.33%) were tri- [71] implemented in the t.test R function. allelic, 1,864,388 (3.18%) had a deletion allele and 57,878 (0.098%) showed four or more alleles (multi-allelic). Measurement of the expression of the ssc‑miR‑326 gene From the set of multi-allelic SNPs, 1765 SNPs mapped to in pigs with different rs319154814 genotypes exons, and 531 of them were clustered in a set of 74 genes Liver samples from nine AA and nine GG pigs enriched in olfactory receptors (see Additional file 7: (rs319154814 genotype) were used to measure the Table S6). This latter result is consistent with the poten- expression of the ssc-miR-326 gene. Total RNA was iso- tial presence of copy number variation (CNV) in such lated with TRIzol as previously described and reverse- locations [73]. Alternative allele frequencies were consist- transcribed with the TaqMan Advanced MicroRNA ently high (> 0.5) for 9.25% of the variants, whereas low Reverse Transcription Kit (Applied Biosystems, Foster (between 0.05 and 0.01) and very low (< 0.01) alternative City, CA), following the instructions of the manufacturer. allele frequencies were detected in 27.85% and 19.62% of Total RNA (10 ng) diluted in a volume of 2 μL was used as SNPs, respectively. a template to generate and pre-amplify cDNA by carry- After filtering, 19,720,314 autosomal whole-genome ing out four consecutive reactions (poly-A tailing, adap- SNPs were selected for assessing population structure tor ligation, reverse transcription and pre-amplification). based on PCA clustering techniques. The PCA showed a The resulting pre-amplified cDNA was diluted to 1:50 in strong genetic differentiation among Asian and European Milli-Q water. To measure ssc-miR-326 expression levels, populations (Fig. 2a). In contrast, domestic pigs and wild both let-7a and miR-26a-5p were chosen as endogenous boars, and particularly those with a European origin, did controls following Timoneda et al. [72]. To conduct the not show such a strong divergence. Asian pigs and wild experiments for the three miRNAs under investiga- boars displayed a certain level of genetic differentiation tion (one target and two controls), available inventoried and were more diverse than their European counterparts. TaqMan advanced miRNA assays (Applied Biosystems, With regard to miRNA variability, the 370 porcine Foster City, CA) for each miRNA were purchased. The miRNA genes annotated in the manually curated miR- reactions were performed in a final volume of 15 μL that Carta v1.1 database [44] were selected, and SNPs within contained: 3.75 μL of 1:50 cDNA, 3 μL of nuclease-free these genes were retrieved. In total, 285 SNPs residing water, 7.5 μL of TaqMan Fast Advanced Master Mix (2×), in 139 miRNAs (37.56% of the total count) were identi- and 0.75 μL TaqMan Advanced miRNA Assay (20×). fied (see Additional file 2: Table S2), which implies that Each reaction was carried out in triplicate in a 384-well the majority of miRNAs are monomorphic. Most of these plate. Reactions were run in a QuantStudio 12K Flex 139 miRNAs (76.98%) included one or two SNPs located Real-Time PCR System (Applied Biosystems, Foster City, within their predicted genomic boundaries, while 18.70% CA). The thermal cycle program was the same as that contained between three up to five SNPs, and 4.32% of used for qPCR experiments on mRNA transcripts. No them displayed more than seven SNPs (Additional file 4: standard dissociation curve was performed. The results Figure S2). Only 43 miRNA SNPs (15.09%) were shared were analyzed applying the ­2−ΔΔCt method [70] for the among the EDM, ADM, EWB and AWB populations relative quantification of miRNA expression (Rq), using (see Additional file 2: Table S2), and showed alternative the group of GG samples as a calibrator. Subsequently, alleles in at least one of the analyzed individuals in each Rq values were ­log2 transformed and the significance of group. The numbers of SNPs segregating in each of the expression differences between GG and AA genotypes four defined groups were 129 (EDM), 201 (ADM), 76 was assessed using a Welch’s t-test for unpaired groups as (EWB) and 172 (AWB), respectively (see Additional file 6: used for mRNA qPCR analyses [71]. Table S5). With regard to precursor and mature regions, 41 and two SNPs were shared among the four populations Results under consideration, respectively (see Additional file 4: Differential segregation of SNPs mapping to microRNA Figure S3A, B). None of the SNPs in the seed regions genes in pigs and wild boars from Europe and Asia were shared by the four porcine populations (see Addi- In total, 58,537,491 million SNPs were identified with the tional file 4: Figure S3C). Only three miRNA SNPs were GATK haplotype caller tool [41] in a dataset comprising found in the European dataset but not in the Asian data- 120 WGS from 40 EDM, 40 ADM, 20 EWB and 20 AWB set. In strong contrast, 55 miRNA SNPs were detected in individuals (see Additional file 1: Table S1). The distri- the Asian dataset but not in the European dataset. Prin- bution of these SNPs within different annotated regions cipal component analyses based on identified autosomal (i.e. protein coding genes, exons, introns, 3′UTRs, miR- miRNA SNPs (N = 260, see Additional file 2: Table S2) NAs, lncRNAs, snRNAs, snoRNAs and pseudogenes) showed poor differentiation between pigs and wild boars Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 9 of 21 Fig. 2 Principal component analysis plots based on SNPs mapping to: a the whole-genome (N = 19,720,314 SNPs), b miRNA genes (N = 285 SNPs), c full 3′UTRs (N = 709,343 SNPs), d 3′UTR 7mer-m8 sites (N = 107,196 SNPs) and e 3′UTR 8mer sites (N = 33,511 SNPs), respectively. The symbol “Q” represents a wild boar sampled from the Middle East (Fig. 2b), whereas the genetic divergence between Euro- precursor regions (see Additional file 2: Table S2), with pean and Asian individuals was still apparent. average alternative allele frequencies of ~ 0.1 for mature When we analyzed population structure based on and seed miRNA polymorphisms. In contrast, the aver- whole-genome autosomal 3′UTR SNPs (N = 709,343 age frequency of the alternative allele observed for SNPs SNPs, Fig. 2c), 3′UTR 7mer-m8 site SNPs (N = 107,196 in precursor areas was ~ 0.15. SNPs, Fig. 2d) and 3′UTR 8mer site SNPs (N = 33,511 Moreover, the observed SNP density adjusted SNPs, Fig. 2e), the genetic differentiation between Asian to 100 bp for miRNA precursor, mature and seed vs. European populations was evident, in close concord- regions consistently followed the order precur- ance with results shown in Fig. 2a, b. However, we also sor > mature > seed when we considered the full set of detected a more pronounced differentiation between 120 WGS. Indeed, ~ one SNP per 100 bp was detected domestic pigs and wild boars, and this observation in precursor regions, whereas ~ 0.86 and ~ 0.42 SNPs was particularly clear for Asian pigs and wild boars per 100 bp were observed in the mature and seed (Fig. 2c–e). regions, respectively (Fig. 4a). After computing a non-parametric comparison of SNP density between Analysis of European and Asian populations shows the mature and seed regions with a Mann Whitney reduced variability in porcine microRNAs U test [50], a significant difference (P-value = 0.016) About 47.76, 57.36, 44.77 and 36.84% of the miRNA SNPs was detected. These results implied a 1.16-fold reduc- showed alternative allele frequencies ≤ 0.1 in the ADM, tion in SNP density between precursor and mature EDM, AWB and EWB populations, respectively (Fig. 3 regions, whereas for seed regions, which are criti- and see Additional file 2: Table S2). Variations located at cal determinants of miRNA–mRNA interactions, mature miRNA and seed regions were enriched in rare or the observed SNP density was ~ 2.4-fold lower than very rare variants compared to the variability of miRNA in precursor regions. This differential distribution Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 10 of 21 Fig. 3 Alternative allele frequency distribution of polymorphisms located at miRNA loci in a Asian domestic pigs (ADM), b European domestic (EDM) pigs, c Asian wild boars (AWB) and d European wild boars (EWB) of the SNP density across miRNA regions (precur- additional anchor pairing to the seed region, showed sor > mature > seed) was also observed in each of the a decreased SNP density, particularly at positions 16 analyses performed in the ADM, EDM, AWB and EWB and 17 (Fig. 4b). In addition, an increased SNP density groups (Fig. 4a). With regard to variants located within was found at positions 19 to 22. Furthermore, the SNP mature miRNAs (N = 64), both inside (N = 12) and density in miRNA genes (N = 370, Fig. 4c) was ~ 2.6- outside seed regions (N = 52), their distribution along fold lower (P-value < 0.01, see “Methods”) than in their the sequence of the mature miRNA (~ 22 nt) showed a flanking regions (± 1 kb). A list of the miRNA SNPs characteristic pattern (Fig. 4b): among all the detected (N = 64) located at mature miRNAs and their genomic SNPs, the first position of the mature miRNA 5′ end, coordinates within the mature miRNA sequence is in which binds to the MID domain of the Argonaute pro- Additional file 8: Table S7. tein in the miRISC complex, exhibited a SNP density of ~ 0.49 SNPs per 100 bp. Such scarcity in polymorphic Statistics of the whole‑genome sequencing of five Duroc sites was also observed when considering the following boars and association analyses positions 2 to 8 in the mature miRNA sequence (seed As previously explained, we sequenced the genomes of region), with an average of ~ 0.42 SNPs per 100 bp five Duroc pigs which founded a population of 350 off- across the whole seed and up to ~ 0.73 SNPs per 100 bp spring with the goal of identifying SNPs in miRNA genes at position 6 of the mature miRNA. In contrast, the and investigating their association with mRNA expres- interval between positions 9 to 12 (with no functional sion and lipid phenotypes. Mean coverage values of the implications in terms of mRNA targeting) showed five pig genomes ranged from 37.67× to 46.6×, with an increased average SNP density of ~ 0.98 SNPs per more than 98.6% of the genome covered by at least 10 100 bp. Positions 13 to 18 of the mature miRNA, which reads in all five samples, and 96.71% of the genome cov- roughly correspond to a functional region providing ered by at least 15 reads. More details about coverage Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 11 of 21 protein coding loci. Fifty-five SNPs and five InDels were located in the 370 miRNA loci annotated in the Sscrofa.11.1 genome assembly according to the miRCarta database [44]. From these, 53 variants (96.36%) were shared with the set of 285 miRNA SNPs found in the 120 WGS reported before (see Additional file 10: Table S9). After performing quality check of the variants mapping to miRNA genes, 49 SNPs and five InDels passed all set filtering criteria. Among the 49 identified miRNA polymorphisms, we selected 12 SNPs on the basis of their location in relevant annotated miRNA loci (Table 1 and see Additional file 10: Table S9) for genotyping and performing association analyses. Microarray expression data comprised 23,998 probes, of which 12,131 (50.55%) and 12,787 (53.28%) corresponded to genes expressed in GM and liver tissues, respectively. To provide a comprehensive view about the association between miRNA genotypes and liver and muscle mRNA levels, we carried out an association analysis between miRNA genotypes and the whole set of expressed probes/genes in these two tissues (Table 2). The results show that there are hundreds of associations at the nomi- nal level, whereas the number of significant associations is strongly reduced after multiple testing correction. This analysis does not take into account the existence of complementarity between the seed of the miRNA and the porcine 3′UTR of target mRNAs or the functional validation data obtained in humans (see “Methods” sec- tion). A second analysis that took these two criteria into account revealed that from the set of 12 genotyped SNPs listed in Table 1, only two SNPs showed significant asso- ciations (q-value < 0.1) with the levels of mRNAs pre- dicted to be miRNA targets (Table 3 and see Additional file 11: Table S10). Further details on the target sites at the 3′UTR of mRNA transcript expression profiles that are significantly associated with miRNA SNP genotypes (Table 3) are shown in Additional file 12: Table S11. All Fig. 4 a SNP density per 100 bp for each analyzed miRNA region considering the full set of 120 porcine whole-genome sequences genes displaying significant associations (q-value < 0.1) (WGS), as well as each WGS from the ADM, EDM, AWB and EWB in this second analysis ranked among the top 1% of sig- groups. b SNP density across mature miRNA regions. 1: anchor nificance (see Additional file 13: Table S12) in the con- (1st 5′ end position), 2: seed (positions 2 to 8) and 3: supplemental dition-free analysis, except for transcripts associated pairing (positions 13 to 18). c SNP density in miRNA loci and flanking with ssc-miR-326 genotypes in the liver, which ranked regions considering a window of ± 1 kb divided in 10 upstream and downstream 100 bp bins among the top 5% (see Additional file 13: Table S12). We observed that for several of the significant associations between miRNA and mRNA pairs reported in Table 3 and see Additional file 11: Table S10, the interaction was and genome mapping parameters are in Additional file 4: classified as 7-mer-m8 by our targeting approach with Figure S4 and Additional file 9: Table S8. After perform- the locate tool of SeqKit toolkit [46] when, in reality, it ing variant calling on mapped reads, 13,839,422 SNPs was of type 8mer with one mismatch (see Additional passed the set quality criteria, whereas 3,721,589 InDels file 12: Table S11). were detected. Moreover, 1,643,861 InDels (44.17%) and 6,034,548 SNPs (43.60%) were located in annotated Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 12 of 21 Table 2 Number of significant associations between the 12 genotyped miRNA SNPs and the whole set of mRNAs expressed in the gluteus medius (GM) skeletal muscle (N = 89) and liver (N = 87) tissue samples from Duroc pigs SNP (chromosome position) Type Significant associations GM LIVER GM LIVER (P-value < 0.05) (P-value < 0.05) (q-value < 0.1) (q-value < 0.1) rs81349391 (1:164,026,014 bp) ssc-miR-339 apical loop (A/G) 456 356 0 0 rs344472188 (2:13,296,736 bp) ssc-miR-130a apical loop (A/G) 320 1252 0 1 rs333787816 (2:65,308,181 bp) ssc-miR-23a precursor stem (T/C) 2075 751 6 1 rs340704946 (4:6,948,743 bp) ssc-miR-30d precursor stem (A/G) 563 398 0 0 rs320008166 (6:56,427,227 bp) ssc-miR-371 precursor stem (T/C) 270 365 2 8 rs323906663 (6:63,491,948 bp) ssc-miR-429 precursor stem (G/A) 328 406 2 8 rs319154814 (9:9,581,989 bp) ssc-miR-326 apical loop (G/A) 372 2400 0 0 rs321151601 (9:39,280,312 bp) ssc-miR-34c mature region (C/A) 947 392 0 0 rs341950320 (12:36,947,510 bp) ssc-miR-378-2 precursor stem (G/A) 792 322 0 0 rs334680106 (13:100,083,229 bp) ssc-miR-15b precursor stem (C/T) 277 205 0 0 rs335924546 (17:10,758,828 bp) ssc-miR-486 precursor stem (C/T) 1232 171 0 0 rs334590580 (18:18,341,582 bp) ssc-miR-335 precursor stem (T/C) 665 219 1 0 This is a condition-free association analysis in which all expressed mRNAs are taken into account regardless of the existence of miRNA–mRNA complementarity and the availability of functional data validating miRNA–mRNA interactions in humans (see “Methods” for more information) Table 3 Top significant associations after correction for multiple testing (q-value < 0.1) between miRNA SNPs and the mRNA expression of their potential targets in the gluteus medius (GM) skeletal muscle (N = 89) and liver (N = 87) tissue samples from Duroc pigs SNP (chromosome Type Tissue Probe ID Gene δa seb P-value q-valuec position) rs333787816 ssc-miR-23a GM Ssc.1790.1.S1_at ENSSSCG00000000019 NUP50 − 0.1052 0.0288 3.849E−04 6.157E−02 (2:65,308,181 bp) precursor Ssc.12493.1.A1_at ENSSSCG00000024027 PAFAH1B2 − 0.2021 0.0565 3.921E−04 6.157E−02 stem (T/C) Ssc.13786.1.A1_at ENSSSCG00000003857 ZYG11B 0.2614 0.0730 4.341E−04 6.157E−02 Ssc.8682.2.A1_at ENSSSCG00000014240 CSNK1G3 − 0.1532 0.0434 4.515E−04 6.157E−02 Ssc.23976.1.S1_at ENSSSCG00000012252 DDX3X − 0.4318 0.1217 5.095E−04 6.157E−02 Ssc.4948.1.S1_at ENSSSCG00000034725 UBE2R2 0.1726 0.0493 5.906E−04 6.157E−02 Ssc.21303.1.S1_at ENSSSCG00000003630 AGO1 0.1016 0.0275 6.775E−04 6.157E−02 Ssc.24035.2.A1_at ENSSSCG00000005935 AGO2 0.2896 0.0864 8.339E−04 6.157E−02 Ssc.3706.1.S2_at ENSSSCG00000040317 SOD2 − 0.2987 0.0896 8.880E−04 6.157E−02 Ssc.3706.1.S1_at ENSSSCG00000040317 SOD2 − 0.4224 0.1261 8.884E−04 6.157E−02 rs319154814 ssc-miR-326 LIVER Ssc.9544.2.S1_a_at ENSSSCG00000016101 CFLAR 0.3246 0.0997 1.264E−03 5.785E−02 (9:9,581,989 bp) apical loop Ssc.11661.2.S1_at ENSSSCG00000009828 PPP1CC 0.4508 0.1436 1.784E−03 5.785E−02 (G/A) Ssc.31172.3.S1_at ENSSSCG00000008601 SDC1 0.2507 0.0808 2.486E−03 5.785E−02 Ssc.11044.1.A1_at ENSSSCG00000003643 SF3A3 0.1808 0.0561 2.598E−03 5.785E−02 Ssc.23242.1.A1_at ENSSSCG00000039426 FSTL1 0.2407 0.078 2.606E−03 5.785E−02 Ssc.12222.1.S1_at ENSSSCG00000013233 CELF1 − 0.1714 0.0574 3.210E−03 5.939E−02 Ssc.10946.2.A1_at ENSSSCG00000011917 NAA50 0.2853 0.0993 4.070E−03 6.454E−02 Ssc.4516.2.S1_at ENSSSCG00000013592 ELAVL1 0.1666 0.0598 6.226E−03 8.639E−02 In this analysis, we consider only the expression of mRNAs with 3′UTR sites complementary to the seeds of the investigated miRNAs and with such interactions having been functionally validated in humans (see “Methods” for more information) a δ: estimated allele substitution effect b se: standard error of the substitution effect c q-value: q-value calculated with the false discovery rate (FDR) approach [65] Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 13 of 21 Two SNPs in the apical loop of ssc‑miR‑326 and in the rs319154814 genotypes). After comparing the mean precursor region of ssc‑miR‑23a are associated expression for each of the three analyzed mRNA tran- with the mRNA expression of several of their putative scripts in both AA and GG genotype groups, only the targets PPP1CC transcripts showed a significantly reduced When we analyzed the association between the expression in AA pigs compared with their GG coun- rs319154814 (G/A) polymorphism located in the apical terparts (P-value = 0.027, Fig. 6a). The other two tran- loop of ssc-miR-326 and gene expression data (Table 3), scripts (CFLAR and SF3A3) also showed a reduced significant results were obtained after multiple testing expression in pigs with genotype AA (Fig. 6a), but correction (q-value < 0.1). More specifically, we detected the results were not significant. We further investi- eight significant associations between this SNP and the gated whether the ssc-miR-326 rs319154814 (G/A) hepatic mRNA expression of experimentally confirmed genotypes were associated with the levels of this very targets of this miRNA (Table 3). For instance, the hepatic same miRNA in the liver by using a specific Taqman mRNA levels of the genes cellular FLICE-like inhibitory Advanced miRNA probe assay. In this qPCR analysis, protein (CFLAR), protein phosphatase 1 catalytic subunit pigs with AA genotypes for the rs319154814 polymor- γ (PPP1CC), syndecan 1 (SDC1), splicing factor 3A subu- phism showed an increased ~ 1.9-fold expression of nit 3 (SF3A3) and follistatin-like 1 (FSTL1) were among ssc-miR-326 transcripts compared with their GG coun- the most significantly associated with ssc-miR-326 terparts, although this difference did not reach statisti- genotypes (Table 3). The expression levels of seven of cal significance (P-value = 0.178, Fig. 6b). these mRNAs associated with ssc-miR-326 rs319154814 The rs333787816 (T/C) polymorphism, located in (G/A) genotypes (Table 3) were reduced in pigs that the precursor region, immediately downstream of the are homozygous for the mutated allele (AA, N = 32), as mature ssc-miR-23a sequence, was also significantly depicted in Fig. 5. associated with 36 experimentally confirmed targeted The mRNA levels of the PPP1CC, CFLAR and SF3A3 genes in the GM muscle tissue (q-value < 0.1, Table 3, mRNA transcripts were further measured by RT-qPCR and see Additional file 11: Table S10). From these, it is in 10 liver samples (5 belonging to each AA or GG worth mentioning the Argonaute RISC component 1 Fig. 5 Hepatic mRNA expression levels of the CFLAR, PPP1CC, SDC1, SF3A3, FSTL1, NAA50 and ELAVL1 genes according to the genotype of the rs319154814 (G/A) apical loop SNP in the ssc-miR-326 gene. The number of individuals representing each genotype were: GG (N = 17), GA (N = 37) and AA (N = 32) Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 14 of 21 Fig. 6 a Barplots depicting qPCR ­log2 transformed relative quantities (Rq) for CFLAR, PPP1CC and SF3A3 mRNA transcripts measured in the liver tissue. The number of individuals representing each genotype were: GG (N = 5) and AA (N = 5). b Barplots depicting qPCR ­log2 transformed relative quantities (Rq) for ssc-miR-326 transcripts measured in the liver tissue. The number of individuals representing each genotype were: GG (N = 9) and AA (N = 9) (AGO1) and the Argonaute RISC catalytic component 2 deacetylase 3 (HDAC3), period 1 (PER1) and the already (AGO2) mRNAs. Both AGO1 and AGO2 genes showed mentioned PPP1CC gene. Moreover, negative regula- lowered mRNA expression levels in homozygous CC tion of DDX58/IFIH1 signaling of type I interferon pigs with respect to their TT and TC counterparts (data response and p53 signaling were pathways enriched in not shown). the sets of mRNAs nominally associated with ssc-miR- 130a genotypes (see Additional file 14: Table S13). Pathway enrichment analyses of significantly associated mRNA targets provide insights of putative metabolic Porcine lipid phenotypes are associated functions of miRNAs with the genotypes of miRNA genes After performing pathway enrichment analyses on the We also evaluated the association between miRNA sets of putative targeted mRNA transcripts that are SNPs and several lipid-related phenotypes recorded in significantly associated with miRNA SNP genotypes the Lipgen population (see Additional file 3: Table S3). (P-value < 0.05), several biological categories were iden- Only the rs319154814 (G/A) variant in the ssc-miR-326 tified and reported in Additional file 14: Table S13. To gene was significantly associated (q-value < 0.1) mention a few of the most relevant findings, the ssc- with lipid traits (Table 4 and see Additional file 15: miR-23a genotype was associated with several regula- Table S14). More specifically, we found significant asso- tory pathways (post-transcriptional silencing by small ciations with myristic acid (C14:0) content in both LD RNAs, competing endogenous RNAs regulate PTEN and GM muscles, as well as with gadoleic acid (C20:1) translation, regulation of PTEN mRNA translation, reg- content and the ratio between PUFA and MUFA in ulation of RUNX1 expression and activity, etc.) involv- the LD muscle (Table 4 and see Additional file 15: ing genes such as AGO1, AGO2 or TNRC6C, which play Table S14). We also observed several additional asso- relevant roles in miRNA-mediated gene regulation [74]. ciations between the rs319154814 (G/A) SNP and fatty Among the sets of mRNAs nominally associated with acid composition traits, but they were significant only the ssc-miR-30d genotype in GM, we identified the reg- at the nominal level (P-value < 0.01, Table 4). Other api- ulation of ornithine decarboxylase activity or cellular cal loop SNPs such as rs81349391 (A/G) at ssc-miR-339 response to hypoxia, as well as mature mRNA transport and rs344472188 (T/C) at ssc-miR-130a were signifi- after splicing (see Additional file 14: Table S13). cantly associated at the nominal level (P-value < 0.01) In the liver tissue, the ssc-miR-326 genotype was with palmitic acid (C16:0) content and SFA and UFA significantly associated with the expression of sev- proportion in the GM muscle, as well as with backfat eral circadian clock regulatory transcripts e.g., histone thickness. As shown in Table 4, a SNP located in the Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 15 of 21 Table 4 Significant associations at the nominal level (P-value < 0.01) and after multiple testing correction (q-value < 0.1; in italics) between 12 genotyped miRNA SNPs and lipid phenotypes recorded in a population of Duroc pigs (N = 345) SNP (chromosome position) Type Trait δ se P-value q-value rs81349391 (1:164,026,014 bp) ssc-miR-339 apical loop (A/G) GM (C16:0) − 0.3234 0.1396 2.005E−02 3.772E−01 GM SFA − 0.4394 0.2030 2.969E−02 3.772E−01 GM UFA 0.4392 0.2030 2.978E−02 3.772E−01 rs344472188 (2:13,296,736 bp) ssc-miR-130a apical loop (T/C) BFTLR − 1.2084 0.5914 4.087E−02 9.683E−01 rs340704946 (4:6,948,743 bp) ssc-miR-30d precursor stem (A/G) GM (C20:0) − 0.0219 0.0098 2.541E−02 5.779E−01 rs320008166 (6:56,427,227 bp) ssc-miR-371 precursor stem (T/C) GM IMF 0.5368 0.2461 3.104E−02 9.060E−01 rs323906663 (6:63,491,948 bp) ssc-miR-429 precursor stem (G/A) GM IMF 0.5682 0.2482 2.348E−02 8.923E−01 rs319154814 (9:9,581,989 bp) ssc-miR-326 apical loop (G/A) LD (C14:0) 0.0885 0.0213 8.722E−05 3.314E−03 GM (C14:0) 0.0700 0.0184 3.942E−04 7.489E−03 LD PUFA/MUFA − 0.0506 0.0188 7.109E−03 7.653E−02 LD (C20:1) 0.0349 0.0125 8.055E−03 7.653E−02 LD (C20:4) − 0.2918 0.1268 2.092E−02 1.438E−01 LD MUFA 0.9512 0.3971 2.485E−02 1.438E−01 LD (C18:1) 0.8349 0.3624 2.814E−02 1.438E−01 LD PUFA − 1.0670 0.5004 3.222E−02 1.438E−01 LD (C18:2) − 0.7266 0.3443 3.405E−02 1.438E−01 GM (C18:0) − 0.2245 0.1115 4.423E−02 1.681E−01 rs321151601 (9:39,280,312 bp) ssc-miR-34c mature region (C/A) LD (C20:0) − 0.0262 0.0119 2.655E−02 8.171E−01 rs335924546 (17:10,758,828 bp) ssc-miR-486 precursor stem (C/T) GM (C20:4) 0.3396 0.1465 2.100E-02 4.555E-01 LD (C14:0) − 0.0586 0.0256 2.937E−02 4.555E−01 LD IMF 0.3300 0.1553 4.748E−02 4.555E−01 rs334590580 (18:18,341,582 bp) ssc-miR-335 precursor stem (T/C) GM (C16:0) − 0.4305 0.1558 5.862E−03 1.201E−01 GM (C20:0) 0.0331 0.0121 6.321E−03 1.201E−01 GM UFA 0.5231 0.2278 2.147E−02 2.045E−01 GM SFA − 0.5229 0.2278 2.153E−02 2.045E−01 GM: gluteus medius skeletal muscle; LD: longissimus dorsi skeletal muscle; BFTLR: backfat thickness at the last rib; IMF: intramuscular fat content; SFA: saturated fatty acids content; UFA: unsaturated fatty acids content; PUFA: polyunsaturated fatty acids content; MUFA: monounsaturated fatty acids content; C14:0: myristic content; C16:0: palmitic content; C18:0: stearic content; C18:1: oleic content; C18:2: linoleic content; C20:0 arachidic content; C20:1: gadoleic content; C20:4: arachidonic content a δ: estimated allele substitution effect b se: standard error of the substitution effect c q-value: q-value calculated with the false discovery rate (FDR) approach [65] precursor 3′ stem of ssc-miR-30d showed a nominally Asian pigs and wild boars are more diverse than their significant association with arachidic acid (C20:0) con- European counterparts and that both gene pools split tent. Other significant associations at the nominal level during the mid-Pleistocene ~ 1.6 to 0.8 Mya. Calabrian were, for instance, those between the rs334590580 glacial intervals probably favored a restricted gene flow (T/C) SNP located at the precursor stem region of ssc- between these two gene pools [75]. The high variability miR-335 and palmitic and arachidic acids content in of Asian populations could be explained by the fact that GM tissue (Table 4). Sus scrofa emerged as a species in Southeast Asia (~ 5.3 to 3.5 Mya) and then dispersed westwards until reach- Discussion ing Europe ~ 0.8 Mya [76]. This initial founder effect, Divergent patterns of variation for microRNA and 3′UTR combined with the occurrence of strong bottlenecks, polymorphisms in Asian and European pigs and wild boars reduced the genetic diversity of European wild boars The PCA revealed the existence of a detectable genetic [75]. differentiation between Asian and European popula- While genetic differentiation between wild boar and tions, with the latter showing reduced levels of diversity pig populations was clearly detected when the whole- compared to the former. Groenen et al. [75] investi- genome SNP dataset was considered (Fig. 2a), this gated the variability of pig genomes and found that was less evident in the PCA based on miRNA SNPs (Fig. 2b), probably because the small number (N = 285 Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 16 of 21 SNPs) of markers used in this analysis limited the reso- this critical functional region are not well tolerated. This lution with which population differentiation could be distribution is similar to the one that we have observed detected. We also found that the degree of population in domestic pigs and wild boars, with 81, 31, 11, 9 and differentiation between Asian domestic pigs and Asian 5 miRNAs harboring one, two, three, four and five SNPs wild boars increased when PCA are built on the basis of (see Additional file 4: Figure S2). Only four and two miR- SNPs located in the whole 3′UTRs (N = 709,343 SNPs), NAs showed a total of seven and ten polymorphisms 3′UTR 7mer-m8 sites (N = 107,196 SNPs) or 3′UTR 8mer within their sequences (see Additional file 2: Table S2). (N = 33,511 SNPs) sites (Fig. 2c–e). Potential effects of We also detected a non-uniform SNP density along 3′UTR SNPs are: modulation of mRNA expression or the sequence of mature miRNAs (Fig. 4b). Gong et al. changes in secondary structure, stability, localization, [80] showed that SNPs tend to concentrate in the mid- translation, and binding to miRNAs and RNA-bind- dle region of the mature miRNA gene rather than at its ing proteins [77]. Purifying selection is less intense in 5′ and 3′ends. Moreover, the same authors described an 3′UTRs than in protein-coding regions, which implies increased SNP density at positions 9 and 15 of the mature that 3′UTRs evolve faster and accumulate a larger frac- miRNA, a result that agrees well with ours (Fig. 4b). tion of recent polymorphisms contributing to population However, we have also identified an elevated number of differentiation [78]. polymorphic sites at positions 11, 19 and 20, a finding that does not match the human data presented by Gong Low SNP density in microRNA genes and lack of uniform et al. [80]. SNP distribution across sites Several of the sites that show a reduced variability in We have found that, in general, miRNA loci have a sub- porcine miRNAs exert critical functions (Fig. 4b and see stantially lower SNP density in their seeds compared Additional file 8: Table S7). For instance, the first nucle- with mature and precursor regions (Fig. 4a). When we otide of mature miRNAs plays an important role in the analyzed the SNP density in miRNA loci and their flank- loading of the mature miRNA on the Argonaute protein ing regions (± 1 kb), a significantly reduced number of to form the miRISC complex [82]. The second to eighth SNPs per 100 bp was found in miRNAs compared with nucleotides in the mature miRNA correspond to the upstream and downstream flanking sequences (Fig. 4c). seed, where we found a consistently reduced SNP den- These results are in agreement with data reported by sity (Fig. 4a, b) compared with other miRNA regions. Omariba et al. [79] and Saunders et al. [5]. This result was expected because this region has a crucial The low variability of miRNA genes, a feature that is role in modulating the interaction between the mature particularly evident in the seed region, is probably due miRNA and its 3′UTR binding sites. Polymorphisms in to the intense effects of purifying selection. Indeed, the the seed region have the potential to disrupt the proper importance of the miRNA seeds is revealed by the high miRNA–mRNA pairing and thus alter biologically rel- conservation of their sequence across species [16, 23, evant regulatory pathways, which tend to be evolution- 48], as this sequence ultimately determines the success arily conserved [48]. Hence, these variants located in the of miRNA–mRNA interactions [2]. In our study, 221, 52 miRNA seeds (see Additional file 8: Table S7) might favor and 12 SNPs were found in precursor, mature and seed the emergence of novel miRNA–mRNA interactions as regions within miRNA loci, respectively (see Additional well as the abolishment of conserved ones, thus modify- file 2: Table S2 and Additional file 4: Figure S3). Gong ing gene regulatory networks. et al. [80] described the existence of 40% polymorphic In contrast with the first eight 5′ positions of the miRNAs in the human genome but only 16% of them miRNA, we found a higher SNP density at positions 9 to displayed more than one SNP. In a more recent study, 12, which do not contribute substantially to miRNA tar- He et al. [81] reported 1879 SNPs in 1226 (43.6%) human get recognition [2] (Fig. 4b). Positions 13 to 16 facilitate miRNA seed regions, and 97.5% of these polymorphisms 3′-compensatory pairing between the mature miRNA had frequencies lower than 5%. These results agree well and targeted 3′UTRs [83], although only at marginal lev- with the overall distribution of the frequency of miRNA els [84]. Nevertheless, in our porcine dataset only posi- SNPs in the European and Asian populations analyzed in tions 16 and 17 showed a reduced SNP density (Fig. 4b). the current work (Fig. 3). He et al. [81] also demonstrated that 1587, 749, 340, 102, 31, and 4 miRNAs harbored zero, one, two, three, four, and five SNPs, respectively, in their seed regions, which reflects that mutations in Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 17 of 21 Polymorphisms in microRNA genes show associations ssc-miR-326 genotype were CFLAR, PPP1CC, SDC1, with the mRNA expression of several of their predicted SF3A3 and FSTL1. The protein encoded by the PPP1CC targets gene belongs to the protein phosphatase PP1 subfam- We investigated the association of 12 genotyped SNPs ily, which is a ubiquitous serine/threonine phosphatase mapping to miRNA genes and hepatic and muscle mRNA involved in regulating multiple cellular processes through expression in Duroc pigs. Several hundreds of significant dephosphorylation signaling. Among these processes, it associations at the nominal level (P-value < 0.05) were is worth mentioning insulin signaling [88], post-trans- detected in both GM muscle and liver tissues when the lational localization of circadian clock components [89] whole set of microarray expression profiles of mRNA and lipids [90, 91] and glycogen metabolism regulation transcripts were considered (Table 2). One the one hand, [92]. Moreover, the CFLAR gene encodes the cFLIP pro- one pitfall of this analysis is that it does not take into tein involved in the inhibition of Fas-mediated apopto- account whether there is molecular evidence supporting sis [93], and the FSTL1 gene plays a role in the immune the existence of direct miRNA–mRNA interactions dis- inflammatory signaling and fibrosis in the liver [94]. playing significant associations. On the other hand, indi- Pigs homozygous for the derived A-allele of the rect interactions between miRNAs and mRNAs might rs319154814 SNP showed a consistent downregulation also exist i.e., a SNP regulating the expression of a miRNA of the hepatic mRNA expression of the CFLAR, PPP1CC, might repress the translation of a given mRNA through SDC1, SF3A3, FSTL1, NAA50 and ELAVL1 genes (Fig. 5), a direct interaction (an event that would not necessarily a result that was confirmed by qPCR for PPP1CC mRNA imply any transcriptional decay, hence being undetect- transcripts (P-value = 0.027, Fig. 6a). Although in the able in our experimental system) and, in turn, this might qPCR experiment the CFLAR and SF3A3 genes also indirectly affect the expression of other multiple mRNAs showed decreased mRNA levels in pigs with homozygous regulated by the translationally repressed mRNA. In spite AA genotypes for the rs319154814 SNP, such associations of this consideration, we decided to restrict our associa- did not reach statistical significance (Fig. 6a). The signifi- tion analyses to miRNA–mRNA pairs with molecular cant downregulation of at least one top putative mRNA evidence of potential interactions because we believe that target (i.e., PPP1CC) might suggest that the rs319154814 with this stringent approach, we are able to minimize the (G/A) variant could increase the repressive activity of chances of detecting spurious false positive results. The ssc-miR-326. Indeed, when we assessed, with a dedicated genomic positions of miRNA-binding sequences in the qPCR assay, the expression level of ssc-miR-326 tran- 3′UTR of targeted mRNAs showing significant associa- scripts in the liver tissue, we found that ssc-miR-326 ­log2 tions (q-value < 0.1, Table 3) with miRNA genotypes are relative expression was increased by ~ 1.9-fold in pigs in Additional file 12: Table S11. homozygous for the A-allele, although this result was not After correction for multiple testing (q-value < 0.1, statistically significant (P-value = 0.178, Fig. 6b). Table 3), only SNPs mapping to the ssc-miR-23a and ssc- As previously discussed, the rs319154814 (G/A) poly- miR-326 genes displayed significant associations with morphism is located in the apical loop region of ssc- the expression levels of some of their putative mRNA miR-326. Although the apical region of a miRNA does targets. More in detail, the rs333787816 (T/C) polymor- not have a function as critical as the seed, polymor- phism in the ssc-miR-23a gene showed significant asso- phisms located in this specific region can have relevant ciations with several putative target mRNAs, among effects on the structural conformation of the pre-miRNA which are AGO1 and AGO2 transcripts, two essential hairpin. Indeed, Fernandez et al. [19] described a muta- components of the miRNA-mediated silencing machin- tion in the apical loop of hsa-miR-30c (G/A) that induces ery [85, 86]. The rs319154814 (G/A) polymorphism in the a steric disruption of the pri-miRNA folding structure apical loop region of ssc-miR-326 also showed a signifi- of the hairpin, hence creating a bulge around the flank- cant association, after multiple testing correction, with ing downstream CNNC motif that facilitates the SRSF3 the hepatic mRNA expression of several of its predicted factor accessibility to the RNA sequence [95]. In other and experimentally-confirmed targets (Table 3). In con- words, SNPs in the apical region can modify the effi- trast, no association with GM muscle mRNA expression ciency with which the Drosha processing machinery was observed. Tissue-specific differences in the expres- is recruited. The rs319154814 (G/A) polymorphism sion of the miRNA or its mRNA targets might explain detected in our study might have structural conse- such an outcome. Indeed, the analysis of the distribu- quences similar to those described for the hsa-miR-30c tion of miRNA expression across human tissues has apical loop variant [19]. shown that only a minority of miRNAs are expressed One clear limitation of our experimental design is that ubiquitously [87]. The hepatic mRNA targets showing we did not test the existence of miRNA–mRNA inter- the most significant association with rs319154814 (G/A) actions and structural miRNA hairpin modifications Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 18 of 21 by using functional assays. Such experiments would expected, this reduced miRNA variability is particularly be necessary to demonstrate our hypothesis that the accentuated in the seed region, a finding that is likely rs319154814 (G/A) polymorphism in the apical loop of explained by the strong effects of purifying selection pre- ssc-miR-326 has a role in promoting the processing and serving the sequence conservation of this critical site. maturation efficiency of the hairpin precursor, thus mod- We detected one SNP in the apical loop of ssc-miR-326 ifying the expression profiles of ssc-miR-326 and its puta- and another one in the precursor region of ssc-miR-23a tive mRNA targets. that are associated with the mRNA expression of several of their putative targets. If confirmed with functional A polymorphism in the apical loop of microRNA 326 assays, these results would reinforce the need of explor- is associated with fatty acid composition traits ing the role of miRNA variation within and outside the In our study, the only miRNA polymorphism showing seed in the fine-tuning of mRNA expression in pigs. significant associations (q-value < 0.1) with lipid-related phenotypes was the rs319154814 (G/A) SNP in the api- Supplementary Information cal loop of ssc-miR-326. This SNP was associated with The online version contains supplementary material available at https://​doi.​ org/​10.​1186/​s12711-​021-​00632-3. myristic fatty acid content in the GM and LD muscles of Duroc pigs, as well as with PUFA/MUFA ratio and gad- Additional file 1: Table S1. List of whole-genome sequences from Euro‑ oleic fatty acid content in LD. To the best of our knowl- pean domestic pigs (EDM, N = 40), Asian domestic pigs (ADM, N = 40), edge, no direct effect of ssc-miR-326 on the metabolism European wild boars (EWB, N = 20) and Asian wild boars (AWB, N = 20). of myristic fatty acid has been described so far, but Additional file 2: Table S2. Single nucleotide polymorphisms located there are reports suggesting that several of the tar- in microRNA genes and their frequencies in European (E) and Asian (A) domestic pigs (DM) and wild boars (WB). gets of this miRNA might be involved in carbohydrate and lipid metabolism pathways [96, 97]. For instance, Additional file 3: Table S3. Means and standard deviations (SD) of fatness and intramuscular fat content and composition traits recorded in the the PPP1CC transcript, one of the predicted targets of gluteus medius (GM) and longissimus dorsi (LD) muscles of 345 Duroc pigs. miR-326, encodes a subunit of protein phosphatase-1, Additional file 4: Figure S1. Sequence homology between 12 porcine which activates acetyl-CoA carboxylase α (ACACA) and miRNAs harboring SNPs that have been genotyped in the Lipgen 6-phosphofructo2-kinase/fructose-2,6-bisphosphatase population (N = 345) and their corresponding human miRNA ortholo‑ gous sequences. Figure S2. Number of SNPs present in each of the (PFKFB) enzymes, the main regulators of fatty acid annotated porcine miRNA loci that display at least one SNP. Figure S3. synthesis and glycolysis, respectively [91]. The protein Venn Diagrams illustrating the sharing of SNPs located at the (a) precursor, phosphatase-1 also activates lipogenic transcription fac- (b) mature and (c) seed regions of miRNA loci among Asian domestic pigs (ADM), Asian wild boars (AWB), European domestic pigs (EDM) and tors such as sterol regulatory element-binding protein European wild boars (EWB). Figure S4. Boxplot distribution depicting 1 (SREBF1), carbohydrate-responsive element-binding whole-genome sequencing statistics for the five Duroc boars that sired protein (MLXIPL) and, moreover, it dephosphorylates the Lipgen population (N = 345). the DNA-dependent protein kinase encoded by the Additional file 5: Table S4. Set of primers used for the qPCR quantifica‑ tion of three mRNAs putatively targeted by ssc-miR-326. PRKDC gene, which is another main determinant of hepatic lipogenesis [91]. In summary, a potential effect Additional file 6: Table S5. Distribution of SNPs mapping to specific genomic annotated regions (i.e., protein coding genes, exons, introns, of the rs319154814 (G/A) SNP on the synthesis or deg- 3′UTR, miRNA, lncRNA, snRNA, snoRNA and pseudogenes) and segregat‑ radation of myristic acid can be envisaged, but such a ing in the set of whole-genome sequences from 120 Asian and European hypothesis still needs to be confirmed at the functional domestic pigs and wild boars. level. Similarly, Cardoso et al. [98] explored the rela- Additional file 7: Table S6. List of porcine genes showing 3 or more multi-allelic SNPs mapping to exonic regions in the set of 120 whole- tionship between miRNA SNPs located in the seed and genome sequences from Asian and European domestic pigs and wild mature regions of bovine miRNA loci and fatty acid boars. composition traits. By combining lipid-related phe- Additional file 8: Table S7. Description of SNPs mapping to porcine notype data with mRNA and protein abundances in a mature miRNA loci (N = 409) and segregating in the set of 120 whole- genome sequences from Asian and European domestic pigs and wild multi-omics approach, they described three miRNA boars. SNPs that were significantly associated with several Additional file 9: Table S8. Whole-genome sequencing statistics for the unsaturated and polyunsaturated fatty acid traits, as five Duroc boars that sired the Lipgen population (N = 345). well as with the expression levels of some of their pre- Additional file 10: Table S9. List of SNPs mapping to 370 porcine miRNA dicted targets [98]. loci and segregating in the 5 Duroc boars that sired the Lipgen population (N = 345). Conclusions Additional file 11: Table S10. Association analyses between SNPs in miRNA genes and the mRNA levels of their predicted targets in the gluteus MicroRNA genes show divergent patterns of variation medius (GM) skeletal muscle and liver tissues of Lipgen pigs. between Asian and European pigs and wild boars and, in general, they display low levels of polymorphism. As Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 19 of 21 Competing interests Additional file 12: Table S11. Genomic coordinates and sequence of the The authors declare that they have no competing interests. seed (miRNA) and its predicted 3′UTR binding site (mRNA) for putative miRNA–mRNA pairs showing significant associations after multiple testing Author details correction (q-value < 0.1). 1 Centre for Research in Agricultural Genomics (CRAG), CSIC‑IRTA‑UAB‑UB, Additional file 13: Table S12. Association analyses between miRNA SNPs Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. 2 Animal Breed‑ and the mRNA expression profiles of the whole set of expressed probes/ ing and Genetics Program, Institute for Research and Technology in Food genes in the porcine gluteus medius (GM) skeletal muscle and liver tissues and Agriculture (IRTA), Torre Marimon, 08140 Caldes de Montbui, Spain. 3 from Lipgen pigs. CNAG‑CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona, Spain. 4 Departament de Ciència Additional file 14: Table S13. Pathway enrichment analyses of the lists of Animal i dels Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra, probes/genes significantly associated at the nominal level (P-value < 0.05) Barcelona, Spain. with miRNA SNP genotypes. Additional file 15: Table S14. Association analyses between miRNA SNPs Received: 16 April 2020 Accepted: 15 April 2021 and fatness and intramuscular fat content and composition traits recorded in the gluteus medius (GM) and longissimus dorsi (LD) skeletal muscles of 345 Duroc pigs. References Authors’ contributions 1. Gebert LFR, MacRae IJ. Regulation of microRNA function in animals. Nat MA conceived this study. MA and RQ designed the experimental protocols. RQ Rev Mol Cell Biol. 2019;20:21–37. coordinated phenotyping recording and contributed to generate microarray 2. Bartel DP. Metazoan microRNAs. Cell. 2018;173:20–51. expression data. EMS performed the DNA extractions and selected the SNPs 3. Chen CYA, Shyu AB. Mechanisms of deadenylation-dependent decay. to be genotyped. EMS performed all bioinformatic and statistical analyses of Wiley Interdiscip Rev RNA. 2011;2:167–83. the data. MGLS and AC conducted qPCR analyses. DG contributed to popula‑ 4. Eichhorn SW, Guo H, McGeary SE, Rodriguez-Mias RA, Shin C, Baek D, et al. tion structure analyses. RT performed whole-genome sequencing of the mRNA destabilization is the dominant effect of mammalian microRNAs five Duroc boars. EMS and MA drafted the manuscript. All authors read and by the time substantial repression ensues. Mol Cell. 2014;56:104–15. approved the final manuscript. 5. Saunders MA, Liang H, Li WH. Human polymorphism at microRNAs and microRNA target sites. Proc Natl Acad Sci USA. 2007;104:3300–5. Funding 6. Mencía A, Modamio-Høybjør S, Redshaw N, Morín M, Mayo-Merino The research presented in the current publication was funded by Grants F, Olavarrieta L, et al. Mutations in the seed region of human miR-96 AGL2013-48742-C2-1-R and AGL2013-48742-C2-2-R awarded by the Spanish are responsible for nonsyndromic progressive hearing loss. Nat Genet. Ministry of Economy and Competitivity. We also acknowledge the support 2009;41:609–13. of the Spanish Ministry of Science and Innovation for the Center of Excel‑ 7. Auyeung VC, Ulitsky I, McGeary SE, Bartel DP. Beyond secondary structure: lence Severo Ochoa 2020–2023 (CEX2019-000902-S) grant awarded to the primary-sequence determinants license pri-miRNA hairpins for process‑ Centre for Research in Agricultural Genomics (CRAG, Bellaterra, Spain). Emilio ing. Cell. 2013;152:844–58. Mármol-Sánchez was funded by a FPU Ph.D. grant from the Spanish Ministry 8. Roden C, Gaillard J, Kanoria S, Rennie W, Barish S, Cheng J, et al. Novel of Education (FPU15/01733). María Gracia Luigi-Sierra was funded with a Ph.D. determinants of mammalian primary microRNA processing revealed fellowship “Formación de Personal Investigador” (BES-C-2017-0024) awarded by systematic evaluation of hairpin-containing transcripts and human by the Spanish Ministry of Economy and Competitivity. Dailu Guan was genetic variation. Genome Res. 2017;27:374–84. funded by a Ph.D. fellowship from the Scholarship Council of China (CSC). The 9. Clop A, Marcq F, Takeda H, Pirottin D, Tordoir X, Bibé B, et al. A mutation authors thank the CERCA Programme of the Generalitat de Catalunya (Barce‑ creating a potential illegitimate microRNA target site in the myostatin lona, Spain) for their support and those who provided publicly available data. gene affects muscularity in sheep. Nat Genet. 2006;38:813–8. 10. Moszyńska A, Gebert M, Collawn JF, Bartoszewski R. SNPs in micro‑ Availability of data and materials RNA target sites and their potential role in human disease. Open Biol. Microarray expression data used in the current study were deposited in the 2017;7:170019. Gene Expression Omnibus (GEO) public repository and are accessible through 11. Liang H, Li WH. Lowly expressed human microRNA genes evolve rapidly. GEO Series Accession Number GSE115484. Phenotypic and genotypic data‑ Mol Biol Evol. 2009;26:1195–8. sets generated and analyzed during the current study are deposited in the 12. Xiong P, Schneider RF, Hulsey CD, Meyer A, Franchini P. Conservation and Figshare public repository available at https://​figsh​are.​com/​proje​cts/​SNPs_​ novelty in the microRNA genomic landscape of hyperdiverse cichlid miRNA/​78690. The whole-genome sequencing dataset from the five Duroc fishes. Sci Rep. 2019;9:13848. boars is available at the Sequence Read Archive (SRA) database (BioProject: 13. Peterson KJ, Dietrich MR, McPeek MA. MicroRNAs and metazoan mac‑ PRJNA626370). Variant Calling Format (VCF) files from the 120 European and roevolution: insights into canalization, complexity, and the Cambrian Asian pigs belonging to domestic breeds and wild boars, as well as from the explosion. BioEssays. 2009;31:736–47. five sequenced Duroc boars are available at the following link: https://​figsh​are.​ 14. Lee CT, Risom T, Strauss WM. Evolutionary conservation of microRNA com/​proje​cts/​VCF_​PIGs/​93140. regulatory circuits: an examination of microRNA gene complexity and conserved microRNA-target interactions through metazoan phylogeny. DNA Cell Biol. 2007;26:209–18. Declarations 15. Luo J, Wang Y, Yuan J, Zhao Z, Jian LU. MicroRNA duplication acceler‑ ates the recruitment of new targets during vertebrate evolution. RNA. Ethics approval and consent to participate 2018;24:787–802. Animal care and management procedures were carried out in accordance 16. Simkin A, Geissler R, McIntyre ABR, Grimson A. Evolutionary dynam‑ with national guidelines for Good Experimental Practices and they were ics of microRNA target sites across vertebrate evolution. PLoS Genet. approved by the Ethical Committee of the Institut de Recerca i Tecnologia 2020;16:e1008285. Agroalimentàries (IRTA). 17. Li J, Zhang Z. MiRNA regulatory variation in human evolution. Trends Genet. 2013;29:116–24. Consent for publication 18. Sun G, Yan J, Noltner K, Feng J, Li H, Sarkis DA, et al. SNPs in human miRNA Not applicable. genes affect biogenesis and function. RNA. 2009;15:1640–51. 19. Fernandez N, Cordiner RA, Young RS, Hug N, Macias S, Cáceres JF. Genetic variation and RNA structure regulate microRNA biogenesis. Nat Commun. 2017;8:15114. Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 20 of 21 20. Jin Y, Lee CGL. Single nucleotide polymorphisms associated with micro‑ 42. O’Rawe J, Jiang T, Sun G, Wu Y, Wang W, Hu J, et al. Low concordance of RNA regulation. Biomolecules. 2013;3:287–302. multiple variant-calling pipelines: practical implications for exome and 21. Larson G, Cucchi T, Dobney K. Genetic aspects of pig domestication. genome sequencing. Genome Med. 2013;5:28. In: Rothschild MF, Ruvinsky A, editors. The genetics of the pig. 2nd ed. 43. Hasan MS, Wu X, Zhang L. Performance evaluation of indel calling tools Wallingford: CABI; 2011. p. 14–37. using real short-read data. Hum Genomics. 2015;9:20. 22. Ramos-Onsins SE, Burgos-Paz W, Manunza A, Amills M. Mining the 44. Backes C, Fehlmann T, Kern F, Kehl T, Lenhof H-P, Meese E, et al. miRCarta: pig genome to investigate the domestication process. Heredity. a central repository for collecting miRNA candidates. Nucleic Acids Res. 2014;113:471–84. 2018;46:D160–7. 23. Penso-Dolfin L, Moxon S, Haerty W, Di Palma F. The evolutionary dynam‑ 45. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA ics of microRNAs in domestic mammals. Sci Rep. 2018;8:17050. sequences to function. Nucleic Acids Res. 2019;47:D155–62. 24. Wei W, Wang G, Cheng Y, Yang R, Song J, Huang S, et al. A miR-511-bind‑ 46. Shen W, Le S, Li Y, Hu F. SeqKit: a cross-platform and ultrafast toolkit for ing site SNP in the 3′UTR of IGF-1 gene is associated with proliferation FASTA/Q file manipulation. PLoS One. 2016;11:e0163962. and apoptosis of PK-15 cells. In Vitro Cell Dev Biol Anim. 2019;55:323–30. 47. Grimson A, Farh KKH, Johnston WK, Garrett-Engele P, Lim LP, Bartel DP. 25. Liu R, Deng D, Liu X, Xiao Y, Huang J, Wang F, et al. A miR-18a binding- MicroRNA targeting specificity in mammals: determinants beyond seed site polymorphism in CDC42 3′UTR affects CDC42 mRNA expression pairing. Mol Cell. 2007;27:91–105. in placentas and is associated with litter size in pigs. Mamm Genome. 48. Friedman RC, Farh KKH, Burge CB, Bartel DP. Most mammalian mRNAs are 2019;30:34–41. conserved targets of microRNAs. Genome Res. 2009;19:92–105. 26. Ma C, Song H, Yu L, Guan K, Hu P, Li Y, et al. MiR-762 promotes porcine 49. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. immature Sertoli cell growth via the ring finger protein 4 (RNF4) gene. Sci Principal components analysis corrects for stratification in genome-wide Rep. 2016;6:32783. association studies. Nat Genet. 2006;38:904–9. 27. Stachowiak M, Szydlowski M, Flisikowski K, Flisikowska T, Bartz M, 50. Mann HB, Whitney DR. On a test of whether one of two random variables Schnieke A, et al. Polymorphism in 3′ untranslated region of the pig is stochastically larger than the other. Ann Math Stat. 1947;18:50–60. PPARA​ gene influences its transcript level and is associated with adipose 51. Vidal O, Noguera JL, Amills M, Varona L, Gil M, Jiménez N, et al. tissue accumulation. J Anim Sci. 2014;92:2363–71. Identification of carcass and meat quality quantitative trait loci in a 28. Bartz M, Szydlowski M, Kociucka B, Salamon S, Jeleń HH, Switonski Landrace pig population selected for growth and leanness. J Anim Sci. M. Transcript abundance of the pig stearoyl-CoA desaturase gene has 2005;83:293–300. no effect on fatty acid composition in muscle and fat tissues, but its 52. Gallardo D, Pena RN, Amills M, Varona L, Ramírez O, Reixach J, polymorphism within the putative microRNA target site is associated et al. Mapping of quantitative trait loci for cholesterol, LDL, HDL, with daily body weight gain and feed conversion ratio. J Anim Sci. and triglyceride serum concentrations in pigs. Physiol Genomics. 2013;91:10–9. 2008;35:199–209. 29. Shao GC, Luo LF, Jiang SW, Deng CY, Xiong YZ, Li FE. A C/T mutation in 53. Gallardo D, Quintanilla R, Varona L, Díaz I, Ramírez O, Pena RN, et al. microRNA target sites in BMP5 gene is potentially associated with fatness Polymorphism of the pig acetyl-coenzyme A carboxylase α gene is asso‑ in pigs. Meat Sci. 2011;87:299–303. ciated with fatty acid composition in a Duroc commercial line. Anim 30. Bartz M, Koscianska E, Szczerbal I, Nowacka-Woszuk J, Kociucka B, Genet. 2009;40:410–7. Salamon S, et al. Polymorphism of the porcine miR-30d is associated 54. Cánovas A, Quintanilla R, Amills M, Pena RN. Muscle transcriptomic with adipose tissue accumulation, its fatty acid profile and the ME1 gene profiles in pigs with divergent phenotypes for fatness traits. BMC expression. Livest Sci. 2015;182:54–7. Genomics. 2010;11:372. 31. Lee JS, Kim JM, Lim KS, Hong JS, Hong KC, Lee YS. Effects of polymor‑ 55. Quintanilla R, Pena RN, Gallardo D, Cánovas A, Ramírez O, Díaz I, et al. phisms in the porcine microRNA MIR206/MIR133B cluster on muscle fiber Porcine intramuscular fat content and composition are regulated and meat quality traits. Anim Genet. 2013;44:101–6. by quantitative trait loci with muscle-specific effects. J Anim Sci. 32. Kim JM, Lim KS, Hong JS, Kang JH, Lee YS, Hong KC. A polymorphism 2011;89:2963–71. in the porcine miR-208b is associated with microRNA biogenesis and 56. González-Prendes R, Mármol-Sánchez E, Quintanilla R, Castelló A, Zidi expressions of SOX-6 and MYH7 with effects on muscle fibre characteris‑ A, Ramayo-Caldas Y, et al. About the existence of common determi‑ tics and meat quality. Anim Genet. 2015;46:73–7. nants of gene expression in the porcine liver and skeletal muscle. BMC 33. Chai J, Chen L, Luo Z, Zhang T, Chen L, Lou P, et al. Spontaneous single Genomics. 2019;20:518. nucleotide polymorphism in porcine microRNA-378 seed region leads to 57. Ramayo-Caldas Y, Mármol-Sánchez E, Ballester M, Sánchez JP, González- functional alteration. Biosci Biotechnol Biochem. 2018;82:1081–9. Prendes R, Amills M, et al. Integrating genome-wide co-association and 34. Choi JW, Chung WH, Lee KT, Cho ES, Lee SW, Choi BH, et al. Whole- gene expression to identify putative regulators and predictors of feed genome resequencing analyses of five pig breeds, including Korean wild efficiency in pigs. Genet Sel Evol. 2019;51:48. and native, and three European origin breeds. DNA Res. 2015;22:259–67. 58. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy—analysis of Affymetrix 35. Kim H, Song KD, Kim HJ, Park W, Kim J, Lee T, et al. Exploring the GeneChip data at the probe level. Bioinformatics. 2004;20:307–15. genetic signature of body size in Yucatan miniature pig. PLoS One. 59. Smedley D, Haider S, Durinck S, Pandini L, Provero P, Allen J, et al. The 2015;10:e0121732. BioMart community portal: an innovative alternative to large, central‑ 36. Frantz LAF, Schraiber JG, Madsen O, Megens HJ, Cagan A, Bosse M, et al. ized data repositories. Nucleic Acids Res. 2015;43:W589–98. Evidence of long-term gene flow and selection during domestication 60. Cayuela JM, Garrido MD, Bañón SJ, Ros JM. Simultaneous HPLC analysis from analyses of Eurasian wild and domestic pig genomes. Nat Genet. of α-tocopherol and cholesterol in fresh pig meat. J Agric Food Chem. 2015;47:1141–8. 2003;51:1120–4. 37. Bosse M, Megens HJ, Madsen O, Crooijmans RPMA, Ryder OA, Austerlitz 61. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, F, et al. Using genome-wide measures of coancestry to maintain diversity et al. PLINK: a tool set for whole-genome association and population- and fitness in endangered and domestic pig populations. Genome Res. based linkage analyses. Am J Hum Genet. 2007;81:559–75. 2015;25:970–81. 62. Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for 38. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for illu‑ association studies. Nat Genet. 2012;44:821–4. mina sequence data. Bioinformatics. 2014;30:2114–20. 63. Karagkouni D, Paraskevopoulou MD, Chatzopoulos S, Vlachos IS, Tast‑ 39. Warr A, Affara N, Aken B, Beiki H, Bickhart DM, Billis K, et al. An improved soglou S, Kanellos I, et al. DIANA-TarBase v8: a decade-long collection pig reference genome sequence to enable pig genetics and genomics of experimentally supported miRNA-gene interactions. Nucleic Acids research. Gigascience. 2020;9:giaa051. Res. 2018;46:D239–45. 40. Li H. Aligning sequence reads, clone sequences and assembly contigs 64. Leclercq M, Diallo AB, Blanchette M. Prediction of human miRNA target with BWA-MEM. 2013. http://​arxiv.​org/​abs/​1303.​3997. genes using computationally reconstructed ancestral mammalian 41. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, sequences. Nucleic Acids Res. 2017;45:556–66. et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303. Mármol‑Sánchez et al. Genet Sel Evol (2021) 53:43 Page 21 of 21 65. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practi‑ 83. Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. cal and powerful approach to multiple testing. J R Stat Soc Ser B. 2009;136:215–33. 1995;57:289–300. 84. Salomon WE, Jolly SM, Moore MJ, Zamore PD, Serebrov V. Single-mole‑ 66. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, cule imaging reveals that Argonaute reshapes the binding properties of et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped its nucleic acid guides. Cell. 2015;162:84–95. gene ontology and pathway annotation networks. Bioinformatics. 85. Ruda VM, Chandwani R, Sehgal A, Bogorad RL, Akinc A, Charisse K, et al. 2009;25:1091–3. The roles of individual mammalian Argonautes in RNA interference 67. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. in vivo. PLoS One. 2014;9:e101749. Cytoscape: a software environment for integrated models of biomo‑ 86. Huang V, Zheng J, Qi Z, Wang J, Place RF, Yu J, et al. Ago1 interacts with lecular interaction networks. Genome Res. 2003;13:2498–504. RNA polymerase II and binds to the promoters of actively transcribed 68. Ballester M, Ramayo-Caldas Y, Revilla M, Corominas J, Castelló A, Estellé genes in human cancer cells. PLoS Genet. 2013;9:e1003821. J, et al. Integration of liver gene co-expression networks and eGWAs 87. Ludwig N, Leidinger P, Becker K, Backes C, Fehlmann T, Pallasch C, et al. analyses highlighted candidate regulators implicated in lipid metabo‑ Distribution of miRNA expression across human tissues. Nucleic Acids lism in pigs. Sci Rep. 2017;7:46539. Res. 2016;44:3865–77. 69. Rio DC, Ares M, Hannon GJ, Nilsen TW. Purification of RNA by SDS 88. Brady MJ, Saltiel AR. The role of protein phosphatase-1 in insulin action. solubilization and phenol extraction. Cold Spring Harb Protoc. Recent Prog Horm Res. 2001;56:157–73. 2010;2010:pdb.prot5438. 89. Schmutz I, Wendt S, Schnell A, Kramer A, Mansuy IM, Albrecht U. Protein 70. Livak KJ, Schmittgen TD. Analysis of relative gene expression data phosphatase 1 (PP1) is a post-translational regulator of the mammalian using real-time quantitative PCR and the 2­ -ΔΔCT method. Methods. circadian clock. PLoS One. 2011;6:e21325. 2001;25:402–8. 90. Liu X, Song B, Li S, Wang N, Yang H. Identification and functional analysis 71. Welch BL. The generalizaton of ‘student’s’ problem when several differ‑ of the risk microRNAs associated with cerebral low-grade glioma progno‑ ent population varlances are involved. Biometrika. 1947;34:28–35. sis. Mol Med Rep. 2017;16:1173–9. 72. Timoneda O, Balcells I, Córdoba S, Castelló A, Sánchez A. Determina‑ 91. Wang Y, Viscarra J, Kim SJ, Sul HS. Transcriptional regulation of hepatic tion of reference microRNAs for relative quantification in porcine lipogenesis. Nat Rev Mol Cell Biol. 2015;16:678–89. tissues. PLoS One. 2012;7:e44413. 92. Toole BJ, Cohen PTW. The skeletal muscle-specific glycogen-targeted 73. Paudel Y, Madsen O, Megens HJ, Frantz LAF, Bosse M, Crooijmans protein phosphatase 1 plays a major role in the regulation of glycogen RPMA, et al. Copy number variation in the speciation of pigs: a possible metabolism by adrenaline in vivo. Cell Signal. 2007;19:1044–55. prominent role for olfactory receptors. BMC Genomics. 2015;16:330. 93. Ram DR, Ilyukha V, Volkova T, Buzdin A, Tai A, Smirnova I, et al. Balance 74. Dueck A, Ziegler C, Eichner A, Berezikov E, Meister G. MicroRNAs associ‑ between short and long isoforms of cFLIP regulates Fas-mediated apop‑ ated with the different human Argonaute proteins. Nucleic Acids Res. tosis in vivo. Proc Natl Acad Sci USA. 2016;113:1606–11. 2012;40:9850–62. 94. Mattiotti A, Prakash S, Barnett P, van den Hoff MJB. Follistatin-like 1 in 75. Groenen MAM, Archibald AL, Uenishi H, Tuggle CK, Takeuchi Y, Rothschild development and human diseases. Cell Mol Life Sci. 2018;75:2339–54. MF, et al. Analyses of pig genomes provide insight into porcine demogra‑ 95. Hargous Y, Hautbergue GM, Tintaru AM, Skrisovska L, Golovanov AP, phy and evolution. Nature. 2012;491:393–8. Stevenin J, et al. Molecular basis of RNA recognition and TAP binding by 76. Frantz LAF, Schraiber JG, Madsen O, Megens HJ, Bosse M, Paudel Y, et al. the SR proteins SRp20 and 9G8. EMBO J. 2006;25:5126–37. Genome sequencing reveals fine scale diversification and reticulation 96. Sebastiani G, Grieco FA, Spagnuolo I, Galleri L, Cataldo D, Dotta history during speciation in Sus. Genome Biol. 2013;14:R107. F. Increased expression of microRNA miR-326 in type 1 diabetic 77. Steri M, Idda ML, Whalen MB, Orrù V. Genetic variants in mRNA untrans‑ patients with ongoing islet autoimmunity. Diabetes Metab Res Rev. lated regions. Wiley Interdiscip Rev RNA. 2020;9:e1474. 2011;27:862–6. 78. Bachtiar M, Jin Y, Wang J, Tan TW, Chong SS, Ban KHK, et al. Architecture of 97. Kefas B, Comeau L, Erdle N, Montgomery E, Amos S, Purow B. Pyruvate population-differentiated polymorphisms in the human genome. PLoS kinase M2 is a target of the tumor-suppressive microRNA-326 and regu‑ One. 2019;14:e0224089. lates the survival of glioma cells. Neuro Oncol. 2010;12:1102–12. 79. Omariba G, Xu F, Wang M, Li K, Zhou Y, Xiao J. Genome-Wide analysis of 98. Cardoso TF, Coutinho LL, Bruscadin JJ, da Silva Diniz WJ, Petrini J, Andrade microRNA-related single nucleotide polymorphisms (SNPs) in mouse BGN, et al. Multi-omics approach reveals mir-SNPs affecting muscle fatty genome. Sci Rep. 2020;10:5789. acids profile in Nelore cattle. Genes. 2021;12:67. 80. Gong J, Tong Y, Zhang HM, Wang K, Hu T, Shan G, et al. Genome-wide identification of SNPs in microRNA genes and the SNP effects on micro‑ RNA target binding and biogenesis. Hum Mutat. 2012;33:254–63. Publisher’s Note 81. He S, Ou H, Zhao C, Zhang J. Clustering pattern and functional effect of Springer Nature remains neutral with regard to jurisdictional claims in pub‑ SNPs in human miRNA seed regions. Int J Genomics. 2018;2018:2456076. lished maps and institutional affiliations. 82. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120:15–20. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions