A total number of 208 Duroc pigs (10 males and 198 females) were used in this study, which were obtained from the whole foundation herd of the Beijing Breeding Swine Center. The pigs were located on the same farm, under similar environmental conditions and an identical standard feeding schedule. All animals were inspected for the presence of open wounds, any illness, or abnormal behavior. All pigs are alive and without genetic modification.

We used the phenol-chloroform method to extract genomic DNA from blood. Genotyping of a total number of 50,703 SNPs across the whole genome was performed using the GeneSeek Porcine 50 K SNP Chip (Neogen, Lincoln, NE, United States ). We performed the following quality control through PLINK (V1.90) software (Purcell et al., 2007) and determined the numbers of SNPs in the following categories: 1) SNPs with minor allele frequencies (MAF) ≥5% and 2) SNPs and individual call rates ≥95%. Only autosomal SNPs, with a total number of 40,070 SNPs, were considered for subsequent analyses. The genotyping module of BeadStudio tool (Illumina, Inc., San Diego, CA, United States) was used to determine the genotypes signal intensity of the individuals, including log R ratio (LRR) and B-allele frequency (BAF).

In the present study, the PennCNV (Wang et al., 2007) algorithm was used to identify porcine CNVs. Based on the hidden Markov model (HMM), this algorithm can detect CNVs from SNPs genotyping data, which includes abundant information including the signal intensity and the population frequency (PFB) at each SNP marker, and the distance between SNPs, based on the Sus scrofa (Sscrofa11.1) genome assembly. To salvage the sample affected by genomic wave, a porcine GC-model file was created by calculating the GC content of the 1 Mb region surrounding each SNP and the -gcmodel option in PennCNV was used for adjustment. After detection of CNVs, PennCNV quality filters were used with the following cutoff values: 1) Standard deviation of LRR < 0.30; 2) BAF drift < 0.01; and 3) Waviness factor value within ± 0.05. To reduce the false positive rate, we acquired a CNV containing three or more consecutive SNPs. Referring to the criteria of (Redon et al., 2006), CNV regions (CNVRs) were determined by aggregating overlapping CNVs identified across all samples, which had to be present in at least two individuals. We divided the CNVRs into three types, including gains, losses, and both types (including gain and loss events).

To verify the CNVs identified by PennCNV, we used the QuantiSNP software (Colella et al., 2007) to analyze the same data set as well. The QuantiSNP algorithm assumes an Objective Bayes hidden-Markov model to improve the accuracy of CNVs identification and mapping, and uses a fixed rate of heterozygosity for every SNP. This CNV calling software was run under default parameters. All CNV calls with a Log Bayes Factor <10 were removed.

Using regression analysis, we assessed the relationship between the numbers of CNVRs and the length of each chromosome. From the results of (Feng et al., 2017), we obtained the segmental duplication (SD) regions of the swine genome, and analyzed the relationship between CNVRs and SD using Chi-squared test. At the same time, we used Chi-squared test to compare gene density between CNV regions and non-CNV regions.

To date, only 20 studies have been focused on genome-wide CNV identification in pigs. Of them, two studies employed Sscrofa9.2, and 17 utilized Sscrofa10.2 genome, respectively. To increase the accuracy of the comparisons among studies, CNVRs located on the Sscrofa9.2 and Sscrofa10.2 assembly were converted into the Sscrofa11.1 genome using NCBI Remap tools.

Swine transcripts and annotations were downloaded from the Ensembl database. According to the position of the CNVRs and genes, we identified the protein-coding genes and miRNA partially or completely overlapping with the CNVRs. The DAVID Bioinformatics Resources was used for function analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG).

Based on the structural relationships between protein-coding genes and CNVRs, we classified the genes into three types, as previously suggested by Woodwark and Bateman (2011). Type I CNV-gene was contained entirely within the CNV. Type II CNV-gene partially overlapped the CNV, which were often disrupted and even with fusion genes formed. Type III genes were those that contained the CNV within the gene. To better understand the biology of the aforementioned three types of CNV genes, we investigated their basic characteristics, selective pressures, and functional annotation. The dN and dS values of the pig/human ortholog were obtained from Ensembl Compare database using PAML. We used the Kolmogorov-Smirnov test to compare dN or dS value among three type genes.

Next, we curated dosage-sensitivity gene list, including the imprinted genes, monoallelically genes. These genes were taken from the Database of Chromosomal Imbalance and Phenotype in Humans using Ensembl Resources (DECIPHER, http://decipher.sanger.ac.uk/index) (Firth et al., 2009), the International Standards for Cytogenomic Arrays (ISCA, http://www.iscaconsortium.org) (Riggs et al., 2012), the Catalogue of Parent of Origin Effects database (Morison et al., 2001; Gimelbrant et al., 2007), and the Geneimprint database (www.geneimprint.com) (Chen and Begcy, 2020). The swine genome contains 21 and 369 imprinted and monoallelically expression genes, correspondingly. Based on data from the Ensembl Genome Compare database, we selected the porcine ortholog gene with human dosage-sensitivity genes. Overall, we established a total number of 1,542 dosage-sensitive genes in swine genome, including 166 imprinted genes, 1,043 monoallelically expressed genes.

The sequenced RNA-seq raw data of 27 adult Duroc tissue types, including retina, pancreas, gut, brain, gall bladder, lung, liver, testes, salivary gland, longissimus dorsi, spinal cord, thyroid, lymph, urinary bladder, spleen, prostate, kidney, adrenal gland, esophagus, stomach, heart, nasopharynx, fat, ovary, breast, placenta, and uterus, were downloaded from NCBI SRA (Sequence Read Archive) database with the BioProject number PRJNA392949 (Zhao et al., 2018). After the QC step conducted using FASTQC, Trimmomatic tools (v3.6), RNA-seq clean data were mapped to the Sus scrofa 11.1 genome release version with Hisat2. To obtain expression levels of all genes in the samples of each of the tissue types, fragments per kilobase of exon model per million mapped reads (FPKM) and counts were calculated using StringTie 1.3.4 and FeatureCounts1.6.0 tools, respectively. We analyzed the difference expression profiles between CNV and non-CNV genes in each tissue. To decrease false positive, we adjusted p-value using the Bonferroni method, which the threshold is 1.85E-3.

We used the tissue specificity index (τ) (Yanai et al., 2005) to grade the scalar measurements of the expression specificity, which ranged from 0 for housekeeping genes to 1 for tissue-specific genes. The index τ is defined as τ = Σ i = 1 N ( 1 − x i ) N − 1   , where N is the number of tissues, and x_i is the expression normalized by the maximal expression value.

We used the miRanda tool (Betel et al., 2008) to predict the target gene regulated by miRNAs. The miRanda algorithm integrated biological knowledge on target rules of mammalian microRNAs. In this study, Tot Score and Tot Energy values set 140 and −20, respectively.

To identify the target-recognition preference of miRNAs overlapped CNVR, we employed a random sampling method, based on the procedure proposed by (Wu et al., 2012). The simulation process included two steps: 1) CNV-miRNAs were randomly selected from all miRNAs in the porcine genome, called pseudo-CNV-miRNAs; 2) Based on the relationship between miRNA and the target genes predicted by miRanda, we marked the relationships target genes and pseudo-CNV-miRNAs or pseudo-non-CNV miRNAs, respectively; 3) For each regulatory type, we re-recorded the number of target genes. Steps (a)–(c) were repeated 10,000 times.

In this study, all statistical analyses, including regression analysis, Kolmogorov-Smirnov test, Wilcoxon rank-sum test, Fisher’s exact test, and Chi-squared test were performed using R4.0.0 software.

We identified a total number of 1,371 CNVs within the autosome genome of Duroc populations (Table 1), whose sizes ranged from 8.37 to 2,838.50 kb. The average and the median sizes were 386.30 and 270.05 kb, respectively. The copy number losses were 28.17 times more frequent than the copy number gains (1,324 losses versus 47 gains). The size of the CN losses and CN gains ranged from 8.37 to 2,838.50 kb and from 33.37 to 998.00 kb, respectively. The median and average sizes of the CN losses (270.05 and 390.39 kb) were slightly larger than those of the CN gains (223.71 and 271.03 kb). The distribution of CNVs size ranges are illustrated in Figure 1. In this Duroc swine population, 205 individuals had CNVs, whereas three individuals were without CNVs, with an average number of CNVs per individual genome of 6.59.

By aggregating the overlapping CNVs, a total number of 211 CNVRs across the autosomes were identified, which covered 118.48 Mb of the swine genome and corresponded to 5.23% of the length of the autosomal sequence. Of these CNVRs, 175 were losses and 32 gains, whereas four contained both events which were within the following ranges 8.34–3,882.48 kb, 19.98–1,035.01 kb, and 596.40–2,516.40 kb, respectively. The CNVRs of losses, gains, and both events had means or medians of 596.30, 277.46, and 1,312.96 or 1,069.49 kb, correspondingly. The distribution of the CNVRs size ranges is depicted in Figure 1. In this study, the loss events were approximately 5.47-fold more common than the gain events.

The numbers of CNVRs in each autosome are presented in Table 1, and the location and characteristics of all CNVRs are displayed in Figure 2. Using regression analysis, we found a significant positive linear relationship between the chromosome sequence length and the number of CNVR located on that (p = 5.13E-4) (Figure 3). Longer chromosomes had higher numbers of CNVRs located on that.

These CNVRs were unevenly distributed across the whole genome. Certain chromosomal regions had a relatively high density of CNVs regions, such as each end of SSC1 and SSC2. We found a positive correlation between the number of CNVs and gene density; 77.82% of the total number of CNV regions were located on protein-coding genes, which were significantly more than non-CNV regions (p = 1.65E-4).

In the genomes of many mammals, SD is a necessary condition and catalyst for CNVs formation. In this study, we found that 21.44% of the CNVR sequences (25.63/119.49 Mb) directly overlapped with porcine SD regions, as obtained from the results of (Feng et al., 2017). Using Chi-squared test, we established significant enrichment of CNVRs in the SD region (p = 3.67E-9). It is noteworthy that CNVs are known to co-occur with SDs, and SDs are more abundant in some locations of the porcine genome. Therefore, based on our findings, we suggest that porcine CNVRs are not uniformly distributed in the genome.

In order to verify the CNVs detected by PennCNV, we utilized the QuantiSNP method to analyze the 50 K SNP data and identify CNVs. After applying the quanlity filtering criteria, we detected 2093 CNVs, and 478 CNVRs were determined by merging overlapping CNVs across all samples. The CNVRs distribution and characteristics on genome was similar to those of PennCNV. To decrease the false disvoery rate or increase the reliability of CNVR, We analyzed the overlap between the CNVRs detected by PennCNV and QuantiSNP. For 211 CNVRs identified by PennCNV, 149 CNVRs also detected by QuantiSNP, account for 70.62%. At the same time, 263 CNVRs obtained by QuantiSNP, also identified by PennCNV, account for 55.02%. The different results might be due to the different algorithms between PennCNV and QuantiSNP.

Of the 211 CNVRs we identified, 30 CNVRs did not include genes as annotated in the ENSEMBL database. The loss CNVRs were overlapping with 1,034 (26 miRNAs), the gain ones with 60 (10 miRNAs), and those involving both types with 41 (3 miRNAs) genes. The total number of genes was 1,135, including 1,096 protein-coding genes (also called CNV-gene) (such as the olfactory receptor gene family, taste receptor gene family, DGAT1, PPARA), and 39 miRNA genes, also called CNV-miRNA, including MIR143, MIR335, and MIRLET7.

Of those 1096 CNV-genes, 145 were dosage-sensitive genes. Using Fisher’s exact test, we found dosage-sensitive gene (p = 3.50E-3) enriched in CNV genes. In addition, 602 CNV-genes (54.93%) were regulated by miRNAs, while 7,872 non-CNV-genes (58.13%) were target genes of miRNAs. Therefore, target genes encompassed more non-CNV than CNV genes (p = 4.18E-2).

Here, we compared the expression profiles of CNV and non-CNV genes in 27 adult Duroc pig tissue types. We found that the expression of CNV genes was significant higher than that of non-CNV genes in the prostate (p = 2.62E-5). Of those CNV genes, 22.99% (252) were tissue-specific (TSI > 0.9). However, 24.94% of the non-CNV genes (3,390) were tissue-specific genes, which were significantly more than CNV genes (p = 4.99E-2).

To evaluate the functional annotation of these CNVRs, we conducted gene ontology (GO) and KEGG enrichment analyses using 1096 CNV-genes. We found 34 GO terms and seven KEGG pathways that were statistically significant (p < 0.05) (Table 2). Of the major GO terms 14 were associated with biological processes, 14 with cellular components, and six with molecular function categories. Significant KEGG pathways included the Hippo signaling pathway (ssc04390), Wnt signaling pathway (ssc04310), taste transduction (ssc04742), and glycerolipid metabolism (ssc00561).

According to definition from Woodwark and Bateman (2011), we identified are 862 type I, 206 type II, and 28 type III CNV genes (Table 3) in this study. On average, type I CNV-genes were shorter (43.86 kb) than type II ones (201.33 kb); type III genes were the longest (461.52 kb). Obviously, type I genes were included within CNVs, whereas type III genes contained CNVs. The Kolmogorov-Smirnov test results showed that the mean or median dN plus dS values of type I genes were higher than those of non-CNV genes (Table 3). Therefore, type I genes tend to be rapidly evolving and to have increased mutation rates. Based on the dN + dS values, we also established that type II and III genes mutated more slowly than non-CNV genes (p = 4.24E-3, p = 5.64E-3, respectively). Additionally, the median or mean dN + dS values of the genes overlapping the loss type CNVRs were higher than those of non-CNV genes (p = 2.2E-16). Our findings evidence that type I genes overlapping loss CNVRs were rapidly evolving. The very low dN/dS values of all CNV genes would show that strong negative selection is acting on them. Those genes were remained due to genetic drift or difficulties to remove on genome.

According to functional annotation, type I genes tend to be involved in the regulation of glucose metabolic processes (GO:0010906), skeletal muscle cell differentiation (GO:0035914), glycerolipid metabolism (ssc00561), Wnt signaling pathway (ssc04310), and taste transduction (ssc04742). On the other hand, type II genes were associated with vesicle-mediated transport (GO:0016192), intracellular protein transport (GO:0006886), and Hippo signaling pathway (ssc04390).

A total number of 39 miRNAs (called CNV-miRNAs) were located in the Duroc CNVRs. The remaining miRNAs were referred to as non-CNV miRNAs. Using Wilcox rank-sum test, we found that the number of CNV-miRNA target genes and binding sites were significantly higher than those of non-CNV-miRNA (p _#TargetGene = 1.92E-2, p _{#BindingSites} = 2.43E-2), respectively. Thus, CNV-miRNAs appeared to regulate more target genes than non-CNV-miRNAs.

To characterize CNV-miRNA target genes, we classified all miRNA target genes into three groups, as described earlier (Wu et al., 2012). The first target genes group had 368 genes regulated exclusively by CNV-miRNAs. Of these, 361 target genes were regulated by one CNV-miRNA, whereas the remaining target genes were regulated by two CNV-miRNAs. The second target gene group included 38 genes regulated by a combination of non-CNV mirRNAs and more than one CNV-miRNA. The third group contained 528 target genes regulated only by non-CNV miRNAs.

To investigate the target-recognition preference of CNV-miRNAs, we used a sampling simulation strategy to identify whether the observed number of target genes for each regulatory type could be expected from random sampling. These simulations provided clues for identifying the regulatory patterns of CNV-miRNAs. We found the number of target genes regulated only by two CNV-miRNAs was significantly higher than the expected after the application of random simulations (p = 3.57E-2). In this study, we found seven target genes that were regulated exclusively by two CNV-miRNAs, such as CRK gene regulated by miR-4331 and miR9817. That is, some genes are preferentially targeted by combination of some CNV-miRNAs. Obviously, the copy number alterations of one miRNA influences that of other miRNAs if their binding sites are co-located in the same UTRs. The dosage of miRNAs should be balanced to synergistically regulate the same genes.

In this study, we performed functional enrichment analyses using three groups of target genes regulated by miRNAs. ErbB signaling pathway (ssc04012) was enriched in the first group of target genes, whereas sphingolipid signaling pathway (ssc04071), NF-kappa B signaling pathway (ssc04064), and Wnt signaling pathway (ssc04310) were significantly enriched in the third group.

Recently, accumulating evidence has indicated the widespread distribution of CNVs in the genome. Furthermore, their involvement in genetic variation, phenotypic diversity, and evolutionary adaptation has been acknowledged as a major contribution (Yim et al., 2010). At least 10% of the human genome is considered to be covered by CNVs (Redon et al., 2006; Wong et al., 2007), and speculations exist that human CNVRs may cover up to 13% of the genome sequences (Stankiewicz and Lupski, 2010). In this study, 211 CNVRs were identified, which accounted for 5.23% of the autosomal sequences. Nevertheless, this figure might be conservative, because some CNVs could not be detected, including small (<10 kb) and large CNVs, which is possibly due to the small sample size and low homology probes, and as well as to limitations of current reference genomes such as sequence gaps.

Here, the abundance of loss CNVR events was approximately 5.47-fold higher than that of gain CNVR events, which is consistent with findings obtained in previous studies on cattle (Hou et al., 2011), goat (Fontanesi et al., 2010) and sheep (Hou et al., 2015). This result might be explained by action of biological factors, as suggested by (Fadista et al., 2010). Initially, non-allelic homologous recombination (NAHR) seemed to be one of the main mechanisms responsible for CNVs formation (Zhang et al., 2009). For example, Locke et al. (2006) suggested losses were under stronger selection than gains. In this respect, (Turner et al., 2008) showed that NAHR tended to generate more loss than gain. Furthermore, our results also confirm that type I CNV-genes overlapping loss CNVRs appear to be rapidly evolving.

We observed that CNVRs tended to have a non-uniform distribution in the porcine genome and were enriched in the gene density and segmental duplications regions. In the human genome, CNVRs were found to be more frequently located in some regions in the genome and chromosomes such as the pericentromeric and the subtelomeric regions (Zarrei et al., 2015). Studies have shown that the non-uniform distribution of CNVRs may arise from nearby repetitive sequences. Moreover, human CNVs were significantly overrepresented in simple tandem repeat sequences (Lupski and Stankiewicz, 2005). In primate genome, CNVs were discovered to occur together with SDs (Kim et al., 2008). Therefore, SDs may promote CNV formation (Dumas et al., 2007; Lee et al., 2008). In addition, human CNVRs were observed to be unusually enriched in protein-coding genes. The elevated gene density of CNVs might have been caused by the retention of duplicated sequences that were of adaptive benefit (Nguyen et al., 2006).

Twenty studies were focused on genome-wide identification of porcine CNVs. SNP genotyping platforms, array-based comparative genomic hybridization (aCGH), and next-generation sequencing were applied in these investigations, and a total number of 16,396 CNVRs were detected, with a total length of CNVR ranging from 9.66 to 560.30 Mb in the different studies. The CNVR distributions established in each of the genome examinations are presented in Table 4. Of the 211 CNVR detected in this research, 191 had been also previously detected in earlier studies. These results indicate that approximately 90% of the CNVRs identified here can be validated by previous investigations, whereas 10% of our findings are original, first detected herein.

The most overlapped CNVR counts (98) were consistent with those reported by (Stafuzza et al., 2019), who detected 3,520 CNVR events based on the SNP chip data of 3,520 Duroc pigs. Additionally, some studies (Chen et al., 2012; Jiang et al., 2014; Wang et al., 2015b; Long et al., 2016; Keel et al., 2019) used SNP chip or NGS platforms to identify CNVRs in Duroc pig populations. A total number of 54, 96, 88, 46, and 75 CNVR were identified in these studies that were overlapping, respectively. These results implied that these overlapping CNVRs contained some Duroc breed genome-specific CNVRs.

Remained reports have lower proportion CNVRs overlapped with our study. The issue of low overlapping rates between different reports was also occurred in CNV studies of other studies. We deem that the following reasons could have contributed to the observed differences. First, the study populations of different breeds have various genetic backgrounds. Many previous studies have also shown the presence of breed/line-specific CNVRs in the genome. (Chen et al., 2012) surveyed CNVs in 18 diverse pig populations and discovered that only 20 CNVRs of the 565 CNVRs were available in more than nine pig populations, whereas most CNVRs (72.9%) were limited to only one pig population. Second, there are differences in the sampling methods and genetic drift events among studies on the same breed. Third, different detected platforms have been used, CGH arrays, SNP genotyping, or NGS. Finally, many structural variations in the genome might have remained undiscovered.

The varying copy number of CNV genes changes gene expression due to altered gene dosage and disruption effects by gene structural variations. If CNV is located in the coding region, it alters the protein function, whereas its location in the regulatory region changes the gene expression level. Dosage sensitivity of the included genes is the most popular hypothesis that attempts to explain pathogenic CNVs. We discovered that dosage-sensitive genes were enriched in the CNV regions in the genome of Duroc pigs. The CNV regions of the human genome may be its most dosage-sensitive regions, in which CNVs are likely to be associated with disease development (Zarrei et al., 2015). However, it is worth emphasizing that changes in gene copy number do not always lead to differences in gene expression. Many factors, such as lack of regulatory elements in duplication event, the chromatin environment, and dosage compensation, might maintain stable mRNA levels.

Structural variation and miRNA are two genetic elements which affect gene expression and regulation. Here, we predicted the potential number of miRNA targets of various genes that were located either in CNV or non-CNV regions. In the genome of a Duroc pig population, we found that miRNAs regulated less CNV-gene than non-CNV gene, but the mean number of miRNA per CNV-gene is similar to that per non-CNV gene. (Felekkis et al., 2011) and Jovelin (2015) demonstrated that the genes located in the CNV regions of the human genome were targeted by more miRNA molecules, and CNV genes had more miRNA-binding sites than non-CNV genes. Similarly to the human genome, miRNA regulates more the CNV-gene in the fruit fly genome than non-CNV gene. However, Jovelin (2015) argued that this principle was not universal. In this previous study, worm and zebrafish showed the opposite pattern and had significantly more miRNAs and target sites per non-CNV genes. Therefore, structure variations such as duplication and deletion do not necessarily lead to increased miRNA target sites for CNV-gene. The distinct results among species could result from functional differences between CNV-gene and non-CNV gene, differential abundance of CNV types, and the accuracy of CNV annotations. The evolutionary interaction between miRNAs and CNVs could have been obscured by interspecies differences.

Previous reports have addressed the impact of CNVs on the phenotypic variation of domestic animals species. For instance, (Clop et al., 2012) supposed that bridging the gap between CNV genotypes and complex phenotypes will be the next genetic challenge. In addition, (Fontanesi et al., 2011) showed that duplication of the ASIP (agouti-signaling protein) locus was associated with a grey coat in the Massese sheep. The majority of CNVRs identified in this study overlapped with pig QTLs. Earlier, (Paudel et al., 2015) hypothesized that copy number variations provided the means for rapid adaptation to different environments during speciation/diversification. Here, we also deem that some genes with CNVs have had a possibly prominent role in the ongoing speciation, and might have impacted certain phenotypes through gene dosage alteration or via a positional effect, in which the structural variant might have altered the genomic landscape of the regulatory elements modulating the expression of these genes.

According to the enrichment analysis results and the already known basic gene function, genes related to some specific biological procession were identified, such as DGAT1, DGAT2, MOGAT2, AGPAT2, FABP1, PPARA, ANGPTL3, NPC2 gene involving fat metabolic (see Table 5). DGAT1 and DGAT2 participate in the regulation of energy synthesis and catabolism, and affect fat metabolism and lipid deposition in tissues. MOGAT2 is critically involved in the uptake of dietary fat by the human small intestine. The roles of FABP1 are related to the lipid metabolism regulation by the PPAR signaling pathway. The NPC2 gene has important functions in the transfer of cholesterol from the human lysosome.

Previous studies on human CNV discovery have reported the presence of copy number variable miRNA genes (Wong et al., 2007; Lin et al., 2008). For example, (Marcinkowska et al., 2011) found that approximately 30% of genome miRNAs were located in the human CNV regions. Additionally, (Ha et al., 2009) discovered that miRNAs had an equilibrating role in genomic dosage phenomena. The results of numerous studies have clearly evidenced the feasibility of using the dysregulation of CNV-ncRNAs as a biological marker for disease screening. In this study, we detected 39 miRNA genes that overlapped with CNVRs, including some miRNAs involved in precursor adipocyte differentiation and lipid deposition, such as MIR143, MIR335, MIR378, and MIRLET7 (Table 5). An earlier study by An et al. (2016) revealed that MIR143 was promoted the adipogenic differentiation of porcine bone marrow-derived mesenchymal stem cells. In another investigation, (Li et al., 2016) evaluated differentially expressed liver miRNAs between Tibetan and Yorkshire pigs and identified differentially expressed miRNAs (MIR335 and MIR378) that participated in the glucose and lipid metabolism. It is noteworthy that (Li et al., 2011) adopted a deep sequencing approach to determine the identity and abundance of miRNAs in swine adipose tissue development and found that MIR143 and MIRLET7 were the miRNAs with the highest expression.

Our present analysis results indicate that in the porcine genome CNV-miRNAs tend to target a higher number of genes than non-CNV-miRNAs with a pattern similar to that in the human genome, earlier established by (Wu et al., 2012). These scientists also found that this regulation model might play important roles in the prevention of CNV-miRNA purification. From an evolutionary viewpoint, certain CNV-miRNAs seem to have beneficial effects on biological processes in organisms. Our further analysis revealed that genes targeted by CNV-miRNAs participate in a wide range of biological responses to environmental factors. Obviously, CNV-miRNAs provide a possibility of increasing regulatory complexity using a strategy that increases the number target genes.