Convergent dwarfism consequences of minipigs under independent artificial selections

Abstract

Background

Currently, diverse minipigs have acquired a common dwarfism phenotype through independent artificial selections. Characterizing the population and genetic diversity in minipigs is important to unveil genetic mechanisms regulating their body sizes and effects of independent artificial selections on those genetic mechanisms. However, full understanding for the genetic mechanisms and phenotypic consequences in minipigs still lag behind.

Results

Here, using whole genome sequencing data of 41 pig breeds, including eight minipigs, we identified a large genomic diversity in a minipig population compared to other pig populations in terms of population structure, demographic signatures, and selective signatures. Selective signatures reveal diverse biological mechanisms related to body size in minipigs. We also found evidence for neural development mechanism as a minipig-specific body size regulator. Interestingly, selection signatures within those mechanisms containing neural development are also highly different among minipig breeds. Despite those large genetic variances, PLAG1, CHM, and ESR1 are candidate key genes regulating body size which experience different differentiation directions in different pig populations.

Conclusions

These findings present large variances of genetic structures, demographic signatures, and selective signatures in the minipig population. They also highlight how different artificial selections with large genomic diversity have shaped the convergent dwarfism.

Background

Pig (Sus scrofa) is one of the domesticated animals important as a model animal in biomedical research. Molecular genetic evidence suggested that pigs emerged in Island Southeast Asia (ISEA) approximately 5.3–3.5 million years ago and spread across Eurasia [1]. Since approximately 10,000 years ago, independent domestication in different regions has led to the diverse pig populations with unique traits. Population-specific characteristics have been independently shaped by artificial selection and selective breeding to meet human needs by leaving traces in genomes.

Minipigs have relatively smaller body sizes than other domestic pig populations. With several advantages, such as easy handling even at full maturity [2] and sharing anatomical similarities with humans [3, 4], they are purpose-bred for biomedical research. Currently, several minipig breeds such as the Bama, Göttingen, Mini-LEWE, Minnesota, Wuzhishan, Yucatan, and Korean (ET and L type) minipigs are available [5,6,7]. In contrast to their common dwarfism features, they have highly different breeding histories [5, 6]. Many minipigs such as Mini-LEWE, Göttingen, Minnesota, and Korean minipigs were developed by outbreeding among different pig breeds. Other breeds such as Bama and Yucatan minipigs were constructed by inbreeding native breeds. They are also derived from different geographical origins. For instance, the Wuzhishan and Bama minipigs are from China [6], while the Yucatan minipig is from the Yucatan peninsula of Mexico [5].

A recent study has revealed that Chinese indigenous minipig breeds with similar geographical distributions show large genetic divergences [6]. This suggests that genetic structures and underlying mechanisms involved in the formation of common features may also be highly divergent among minipig breeds due to different breeding histories and geographical origins. However, a full understanding of genetic diversity and its phenotypic consequences in minipigs still lags behind. The main reason is the limited number of pig breeds and target genomic regions used in previous studies. Although genomic diversity among breeds can be better understood when a larger number of breeds are investigated together, many studies have been performed with only one or two target minipig breeds [7,8,9,10]. A recent study [6] has analyzed genomes of five different breeds. However, it was limited in that it used only repeat markers, which are not enough to fully understand their genomes [11]. Hence, many genes such as AR [8], LCORL [12], NR6A1 [12] and VRTN [13] having potential effects on body size in pigs have been revealed. However, whether their effects are universal or breed-specific remain unclear.

Thus, we performed a comprehensive population-level genomic analysis for eight different minipig breeds (three inbred and five outbred breeds, n = 52) with 33 other pig breeds (n = 158) and five outgroup species (n = 6). These populations covered three different domestic pigs (Minipigs, Asian pigs, and European pigs) and two geographically different wild boars. Using whole genome sequencing data of these large populations, distinct genomic architectures and selection signatures of the minipig populations, especially related to body size, were identified.

Methods

Variant calling, evaluation, and annotation

For variant calling, we first collected whole genome sequencing data of 216 samples of 41 pig breeds and five outgroup Sus species (Supplementary Table 1). The quality of collected whole genome sequencing data was examined using FastQC (v.0.11.9) [14]. Low-quality reads and adaptor sequences in the reads were removed using NGStoolkit IlluQC.pl (v.2.3.3) [15]. When there were no adapter sequences in the reads, “N A -s 20 -l 70” parameters were used for NGStoolkit to filter only low-quality reads with low-quality base percentages (Phred quality score < 20) larger than 30%. If adapter sequences were present, “2 A -s 20 -l 70” parameters were used. Low-quality bases (Phred quality score < 20) were next trimmed at the 3' end of the reads. Trimmed reads shorter than 45 bases were then removed using NGStoolkit TrimmingReads.pl with “-q 20 -n 45” parameters. If adapter sequences still remained when checked with FastQC, additional trimming was performed using TrimGalore (v.0.6.0) [16] with the following parameters: “-q 20 --length 45 --paired --illumina”.

Cleaned reads for each sample were mapped against the pig reference assembly (Sscrofa11.1) using BWA MEM (v.0.7.17) [17]. Duplicated reads were marked by Picard MarkDuplicates (v.2.17.11; http://broadinstitute.github.io/picard). Base quality of the reads was then recalibrated using GATK (v.4.1.7.0) [18] BaseRecalibrator and ApplyBQSR. GATK HaplotyeCaller was used to call variants for each sample. The called variants of each sample were combined into a GVCF file using GATK CombineGVCFs and jointly genotyped with GATK GenotypeGVCFs using “--dbsnp” parameter with single-nucleotide polymorphism (SNP) data of pig obtained from the dbSNP database (dbSNP build 150). High-quality SNPs were next collected using GATK SelectVariants with hard-filtering criteria (QD < 2.0 || QUAL < 30.0 || SOR > 3.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0) and VCFtools (v.0.1.13) [19] with the following criteria: “--maf 0.01 --max-missing 0.1”. Finally, SNPs on autosomes and X chromosome were filtered for use in further analyses. To evaluate the quality of the remaining high-quality SNPs, transition-to-transversion ratio (Ts/Tv) was calculated for each breed using VCFtools.

Between-breed analyses

Principal component analysis (PCA), phylogenetic inference, and admixture analysis were performed using filtered SNP data. PCA was performed using the Genome-wide Complex Trait Analysis (GCTA) tool (v.1.91.4) [20]. The genetic relationship matrix for pairs of individuals was estimated using “--make-grm --autosome --autosome-num 18” parameters. Based on the matrix, eigenvalues and eigenvectors were calculated using the “--pca 3” parameter in GCTA. The results were visualized using the ggplot2 R package [21]. A phylogenetic tree was constructed with the maximum likelihood method using SNPhylo [22] with “-A -b” parameters and 100 bootstrap samples. The ancestry of pig individuals was inferred by ADMIXTURE (v.1.3.0) [23] with 200 bootstrap replicates and the number of ancestral clusters K ranged from 2 to 5. The inferred ancestry for each cluster count was visualized using Clumpak [24]. To calculate pairwise fixation index (Fst) scores between breeds, we merged SNPs of two breeds into a GVCF file using GATK CombineGVCFs. Fst scores of each variant in the combined GVCF file were next calculated using VCFtools (v.0.1.13), and they were averaged to obtain pairwise Fst score between two breeds. In this analysis, only breeds with at least three individuals were used to avoid possible biases resulting from small sample sizes. Thus, the Wuzhishan minipig was excluded.

Admixture analyses

We calculated $f$ and $D$ statistics using ADMIXTOOLS (v.7.0.2) [25] with default parameters. The ${f}_{3}$ and $D$ statistics were calculated to test for evidence that each minipig breed was derived by admixture between Asian wild boar (AWB) and European wild boar (EWB). The ${f}_{3}$ statistics were calculated by ${f}_{3}(X;AWB, EWB)$. In the case of D statistics, two different configurations for gene flow from different wild boar populations, $D(EWB, X;AWB, OG)$ and $D(AWB, X;EWB, OG)$, were used. Here, $OG$ indicates an outgroup species (Sumatran wild boar in this analysis). The ${f}_{4}$ ratio statistics were used to quantify the ancestral proportion of an admixture event. For the analysis, the EWB population was randomly divided into two different sub-groups, $EWBa$ and $EWBb$. The ancestral proportion was then estimated by ${f}_{4}\left(EWBa, OG;X,AWB\right) / {f}_{4}\left(EWBa, OG;EWBb,AWB\right)$. Sumatran wild boar was also used as an outgroup species in this analysis. In the above measures, $X$ indicates the target minipig breed.

Within-breed analyses

To confirm demographic signatures within the minipig population, we used runs of homozygosity (ROH), nucleotide diversity ($\pi$), and linkage disequilibrium decay (LD decay). The ROH for each pig individual was calculated using Plink (v.1.90) [26]. In this step, missing genotyped SNPs were imputed using Beagle (v.5.4) [27]. The ROH with a minimum length of 1 Mbp were then obtained with the following parameters used in previous studies [28, 29]: “--homozyg-density 50, --homozyg-gap 1000, --homozyg-kb 1000, --homozyg-snp 100, --homozyg-window-het 5 --homozyg-window-missing 5 --homozyg-window-snp 50 --homozyg-window-threshold 0.05”. The LD decay and π were calculated using SNP data of each breed. Breeds with at least three individuals were used for reliability as in between-breed analyses. Thus, the Wuzhishan minipig was excluded in the analyses. The $\pi$ was estimated for each breed using VCFtools (v.0.1.17) with a 10 Kbp sliding window and a 5 Kbp step size. To estimate LD decay along with varying genomic distances in each breed, squared correlation coefficients (r²) for SNP pairs within 5 Mbp were calculated using PopLDdecay (v3.41) [30] with the following parameters: “-MaxDist 5000”. The rate of LD decay was presented as the genomic distance where r² dropped to half of its maximum value. LD decay patterns and π distributions were visualized using the ggplot2 R package [21]. For additional LD examination, 100 K SNPs were randomly selected across all chromosomes. Pairwise LDs among SNPs in the same chromosome were calculated using plink (v.1.90).

Identification of selective sweeps during domestication

To detect genome-wide selection signatures between different minipig sub-populations, genome-wide ZFst and $\pi$ ratio were calculated with a window size of 10 Kbp and a step size of 5 Kbp. The genome-wide ZFst for each sub-population was calculated against the wild boar (WB) population using VCFtools (v.0.1.17). To calculate the $\pi$ ratio for each window, $\pi$ was calculated for WB (${\pi }_{WB}$), MP1 (${\pi }_{MP1}$), and MP2 (${\pi }_{MP2}$) population using VCFtools (v.0.1.17). The $\pi$ ratio for each window was then calculated by $\text{log}({\pi }_{MP1}/{\pi }_{WB})$ and $\text{log}({\pi }_{MP2}/{\pi }_{WB})$ for MP1 and MP2 sub-populations, respectively.

The sliding genomic regions with high ZFst values (> 2.5) and the bottom 5% $\pi$ ratio in each sub-population (< -0.333 for MP1 and < -1.523 for MP2) were considered as candidate regions under selection. To alleviate geographical biases, we excluded regions with high ZFst values (> 2.5) between AWB and EWB or those in the X-sweep region [8].

To evaluate candidate selective sweep regions, Tajima’s D scores were further calculated with a window size of 10 Kbp and a step size of 5 Kbp for each minipig sub-population and wild boar population. Tajima’s D scores were calculated using VCFtools with “--TajimaD 10000” parameters for non-overlapping 10 Kbp genomic windows.

Genes that overlapped with selective sweep regions were considered as candidate genes under selection. Functional enrichment analyses for candidate genes were performed by g:Profiler [31] with a g:SCS threshold of 0.05.

Differentiation analysis of genes regulating body size

To detect selective sweep regions related to body size regulation in pigs, genome-wide ZFst for minipig sub-populations was first calculated against EDP population using VCFtools (v.0.1.17). To determine the differentiation direction of genes regulating body size of selective sweeps, we calculated $\pi$ ratio, XP-EHH, Tajima’s D scores, and heterozygosity. To calculate $\pi$ ratio for the whole genome, $\pi$ was calculated for EDP (${\pi }_{EDP}$), MP1 (${\pi }_{MP1}$), and MP2 (${\pi }_{MP2}$) populations using VCFtools with a 10 Kbp sliding window and a 5 Kbp step size. The $\pi$ ratio for each window was then calculated by $\text{log}({\pi }_{MP1}/{\pi }_{EDP})$ and $\text{log}({\pi }_{MP2}/{\pi }_{EDP})$ for MP1 and MP2 sub-populations, respectively. We considered genes in the sliding window at the bottom 5% $\pi$ ratio in both MP1 and MP2 sub-populations (< -0.170 for MP1 and < -1.180 for MP2) as MP-specific genes and at the top 5% $\pi$ ratio in both MP1 and MP2 sub-populations (> 0.778 for MP1 and > 0.370 for MP2) as EDP specific genes. Other genes were defined as genes that affect both minipig and EDP populations. XP-EHH scores for MP1 and MP2 sub-populations against the EDP population were estimated and normalized using selscan (v.2.0.0) [32] with default parameters. Normalized XP-EHH scores were averaged in a 10 Kbp sliding window with a 5 Kbp step size. Tajima’s D scores were calculated using VCFtools with “--TajimaD 10000” parameters for non-overlapping 10 Kbp genomic bins for populations. Heterozygosity was defined as $2p\left(1-p\right),$ where $p$ indicates a reference allele frequency, and it was calculated for a 10 Kbp sliding window with a 5 Kbp step size using an in-house Perl script.

Transcriptome analysis

For RNA sequencing data of three tissues (brain, muscle, and liver) for Korean and Bama minipigs, Duroc, and Landrace, Trimmomatic (v.0.36) [33] and SortMeRNA (v.4.2.0) [34] were used with default parameters to remove low-quality reads and rRNA sequences. Filtered RNA-seq reads were mapped to the pig reference genome (Sscrofa11.1) using RSEM (v.1.3.0) [35] with the “--star” parameter and pig reference gene annotation (release 100). RSEM was also used to calculate gene expression levels of each pig breed. Using gene expression values, differentially expressed gene analysis was conducted for all breed pairs for each tissue using the DESeq2 R package (v.1.22.2) [36]. In this analysis, batch factors were added to the model used in DESeq2 to balance samples across experimental batches and control batch differences. Genes with |log2 fold change|≥ 1 and adjusted p-value < 0.05 were identified as differentially expressed genes.

Results

SNP identification

We obtained whole genome sequencing data of 216 samples of eight minipigs, 21 Asian and European domestic pigs (ADP and EDP respectively), 12 Asian and European wild boars (AWB and EWB respectively), and five outgroups (Fig. 1a, Supplementary Table 1; Methods). An average of 15,489,443 single nucleotide polymorphisms (SNPs) were identified in each breed (Supplementary Table 2). The European population (EDP and EWB) had a relatively smaller number of SNPs (11,494,210 SNPs on average) than others (18,300,902 SNPs on average). In all breeds except outgroups, more than 80% of SNPs were found in the dbSNP database (build 150). The transition-to-transversion ratio (Ts/Tv) of SNPs ranged from 2.43 to 2.68, which was comparable to that of a previous study [37].

Fig. 1

Geographical distribution and population structure of pig breeds. a Geographical distribution of pig breeds and outgroups used in this study. The background map was drawn using ggplot2 [21]. Different colored dots indicate different pig groups (Red: European domestic pig; Pink: European wild boar; Blue: Asian domestic pig; Sky: Asian wild boar; Purple: minipig, Grey: Outgroup). b PCA results of 216 samples of pig breed samples and outgroups by the second and third principal components. MP1 and MP2 indicate two sub-populations of minipigs. c Estimated phylogenetic relationships among pig populations and outgroups. In b and c, different colors represent different pig breeds and outgroups as shown in the legend. d Results of admixture analysis with K = 4. Top boxes with different colors indicate different pig populations and outgroup. e D statistics obtained by comparing to Asian and European wild boar (AWB and EWB, respectively). Green and yellow colors represent MP1 and MP2 sub-populations, respectively. Circle and triangle marks show inbred and outbred breeds, respectively. Negative values on the left and right panel indicate significant genetic exchanges of the minipig with AWB and EWB, respectively. f Proportions of Asian and European ancestries estimated using the ${f}_{4}$ ratio statistics. Red and blue colors indicate European and Asian ancestries, respectively

Genetic structure of the minipig population

Using the identified SNPs, the genetic structure of the minipig population was characterized by several analyses such as principal component analysis (PCA), phylogenetic analysis, and admixture analysis (Methods). In PCA, the first three principal components (PCs) explained about 35% of the total genetic variation, and a clear separation of outgroups from others was observed by the first PC (Supplementary Fig. 1). When comparing the second and third PCs (Fig. 1b), commercial pigs (ADP and EDP) and wild boars (AWB and EWB) were largely clustered into two groups, Asian (ADP and AWB) and European (EDP and EWB) pig populations, as reported in a previous study [1]. Interestingly, the minipig population showed a large genetic variance compared to other populations. For example, by the second PC, the minipig population was separated into two sub-groups, one for Bama, Göttingen, Mini-LEWE, Wuzhishan, and Korean minipigs (hereafter called MP1 sub-population), and the other for Minnesota and Yucatan minipigs (hereafter called MP2 sub-population). MP1 and MP2 sub-populations were very close to Asian and European pigs, respectively, based on the second PC. By the third PC, the MP1 sub-population was further distinguished from Asian pigs, and Bama minipigs were separated from other minipigs in the MP1 sub-population. Similar patterns were also observed in estimated phylogenetic relationships and pairwise fixation index (Fst) scores (Fig. 1c, Supplementary Fig. 2).

The admixture analysis showed consistent patterns with larger variability among breeds even in the MP2 sub-population based on varying ancestral cluster number, K (Fig. 1d, Supplementary Fig. 3). Representatively, at K = 4 (Fig. 1d), the MP1 sub-population mainly shared identical ancestry with Asian pigs, whereas the MP2 sub-population was more similar to European pigs. The variability among breeds was also larger in the MP2 sub-population. For example, while Yucatan minipigs had only an EWB ancestor, the majority of Minnesota pigs had an EDP ancestor.

Admixture patterns in minipig populations were further investigated using ${f}_{3}$, $D$, and ${f}_{4}$ ratio statistics [25] by comparing allele-sharing patterns among pig breeds (Methods). Although ${f}_{3}$ scores of all outbred minipigs except Minnesota minipigs were lower than those of inbred minipigs, there was no clear evidence of admixture between Asian and European populations in any minipig breeds because all minipig breeds had positive ${f}_{3}$ scores (Supplementary Fig. 4, Supplementary Table 3). When D statistics were used, MP1 minipigs showed a clearer possibility of genetic flow from the Asian ancestry than MP2 minipigs (larger negative MP1 values; the left panel of Fig. 1e, Supplementary Table 3) and MP2 minipigs showed more obvious evidence of genetic flow from the European ancestry than MP1 minipigs (larger negative MP2 values; the right panel of Fig. 1e, Supplementary Table 3). Consistent patterns were also observed based on ${f}_{4}$ ratio statistics (Fig. 1f, Supplementary Table 3). Furthermore, the MP1 inbred did not show a significant mixture of Asian and European ancestry, whereas MP1 outbred minipigs carried an admixed pattern with a different proportion of the European ancestry ranging from 9.69% (Göttingen) to 18.73% (Mini-LEWE). These results indicate that the minipig population could be largely divided into two different sub-populations carrying different genomic admixtures of Asian and European ancestry.

Demographic signatures of the minipig population

The diversity of demographic signatures within minipigs was examined next using runs of homozygosity (ROH), nucleotide diversity ($\pi$), and linkage disequilibrium decay (LD decay; Methods). A total of 45,875 ROH were identified in 210 individuals of pig breeds except for outgroups (Supplementary Table 4). The majority of the ROH were very small (1–2 Mbp) in all populations, but the frequency of ROH sizes was highly variable within and among populations (Fig. 2a). For Asian and European pigs, ROH patterns consistent with previous studies [38, 39] were identified, such as larger ROH in domestic pigs than their counterpart wild boars and larger ROH in European pigs than in Asian pigs. ROH size distribution was different among minipigs and compared to other pigs. The ROH of crossbred breeds was mainly composed of shorter ROH (1–4 Mbp), and especially, most ROH in Korean minipigs were substantially short (83.60 $\pm$ 2.95% of 1–2 Mbp ROH). Meanwhile, inbred minipigs showed an enrichment of long ROH ($\ge$ 4 Mbp). Especially, the proportion of ROH longer than 8 Mbp was the largest in Bama minipigs (Total length: 976.18 Mbp), and the longest ROH (41.46 Mbp) was also found in the Bama minipig (Supplementary Table 4). Interestingly, although the Minnesota minipig is crossbred, it contained a relatively large proportion of longer ROH, especially with lengths of 2–8 Mbp. The ROH length can be used to infer an inbreeding history. For instance, short and long ROH indicate ancestral and recent inbreeding, respectively. Therefore, the variability of ROH in minipigs suggests that different timing and types of breeding have created discernable demographic signatures in them.

Fig. 2

Different demographic histories of the eight minipig breeds. a Proportion of total sum of ROH (SROH) in different ROH size ranges. Vertical bars represent standard deviations of proportions. ADP: Asian domestic pigs; EDP: European domestic pigs; AWB: Asian wild boar; EWB: European wild boar. b Distribution of SROH and the number of ROH (NROH). c Distribution of nucleotide diversity ($\pi$). d Linkage Disequilibrium decay (LD decay) represented by mean correlation (r²) between SNP pairs in different physical distances

Minipigs were further compared based on the number of ROH (NROH) and the sum of the length of ROH (SROH; Fig. 2b). Here, we assumed that Asian and European pigs are non-admixed populations. Consistent with a previous study [40], outbred minipigs except Minnesota minipigs had fewer ROH and lower sums than inbred minipigs. However, even those outbred minipigs showed larger NROH and SROH than most of the Asian pigs and some European pigs, although the relatively less non-admixed population tends to have larger NROH and SROH than admixed ones [40] (Supplementary Fig. 5). These findings indicate that most of the minipig breeds have undergone population bottleneck which tends to increase ROH size. Especially, Bama and Minnesota minipigs had much larger NROH and SROH than other minipigs, suggesting a stronger population bottleneck in those breeds.

Similar patterns were also obtained in nucleotide diversity ($\pi$) and linkage disequilibrium decay (LD decay) analysis (Fig. 2c, d; Methods). Extensive genome-wide LD and low $\pi$ were consistently observed in Bama and Minnesota minipigs, supporting their strong population bottleneck. Crossbred breeds in MP1 minipigs having a large proportion of Asian ancestry showed similar $\pi$ to the ADP population. The slower LD decay of outbred MP1 breeds than the ADP population also supports population bottleneck in those breeds (Supplementary Table 5). Moreover, the LD decay rate was highly different among outbred MP1 minipigs. Especially, Mini-LEWE minipigs showed a very low LD decay rate similar to Bama and Minnesota minipigs. In the case of the MP2 sub-population, the $\pi$ value was similar, but the LD decay rate was highly different within the sub-population. For an additional examination, we randomly selected 100 K SNPs across all chromosomes and calculated pairwise LDs among SNPs in the same chromosome (Supplementary Fig. 6). Interestingly, the LD among intrachromosomal regions was very highly different among minipig breeds. Especially, Mini-LEWE minipigs retained LDs larger than 0.5 across whole chromosomes. We also identified long-range linkage disequilibrium (LRLD), which is the LD between regions widely separated on a chromosome, in several Mini-LEWE chromosomes.

Selection signatures by domestication in the minipig population

Previous studies [41, 42] have revealed selective sweeps caused by artificial selection during domestication by comparing domestic pig breeds to their wild counterparts. A similar analysis was performed to identify genomic regions under artificial selection in different minipig populations based on Z-transformed fixation index (ZFst) and $\pi$ ratio calculated against AWB and EWB populations (collectively called WB population). To characterize origins of selection signatures, ZFst values for ADP and EDP populations against the WB population were also compared (Methods).

A total of 327 and 1,206 candidate regions under selection (ZFst > 2.5 and $\pi$ ratio < bottom 5%) containing 51 and 343 genes were identified for MP1 and MP2 sub-populations against the WB population, respectively (Fig. 3a and Supplementary Table 6). Candidate regions under selection showed larger absolute differences of Tajima’s D scores against WB population in the corresponding minipig sub-population compared to other sub-populations, providing evidence of sub-population-specific selective sweep signals on candidate regions (Supplementary Fig. 7).

Fig. 3

Results of genome-wide selective sweep analyses for minipig sub-populations against the wild boar (WB) population. a Plots of Z-transformed fixation index (ZFst) and $\pi$ ratio against the WB population for MP1 (top panel) and MP2 sub-populations (bottom panel) in a 10 Kbp sliding window with a 5 Kbp step size. Red dots denote candidate selective sweep regions with ZFst cutoff at 2.5 and $\pi$ ratio cutoff at bottom 5%. b Intersections of selective sweep regions against the WB population for minipig sub-populations and commercial pig populations. ADP: Asian domestic pigs; EDP: European domestic pigs. c ZFst values of selective sweep regions commonly identified with Asian and European domestic pig populations in MP1 (green dots) and MP2 sub-populations (yellow dots). d ZFst distributions of different populations against the WB population in MP1 and MP2 sub-population-specific selective sweep regions

Consistent with results of admixture analyses between Asian and European ancestry, MP1 and MP2 sub-populations shared more selected regions with ADP and EDP populations than other populations, respectively (MP1: 267 vs. 6, MP2: 0 vs. 238; Fig. 3b). ZFst values for most of the common candidate regions in the MP1 or MP2 sub-population were not significantly higher than those of their corresponding pig population (Fig. 3c), indicating that alleles in those selective sweep regions mainly came from the corresponding pig (ADP or EDP) population. MP1- and MP2-specific selective sweep regions (54 and 967 total, respectively, as shown in Fig. 3b) did not show significant ZFst values when tested in other sub-populations against the WB population (Fig. 3d), suggesting that they were mainly shaped by minipig sub-population-specific selection pressures.

Gene ontology (GO) enrichment analysis identified several biological process (BP) GO terms related to actin protein for the 343 genes in the MP2-specific selective sweep regions (Supplementary Table 7). Also, in cellular component (CC) GO terms, genes related to cell or neuron projection associated with the actin protein [43] (Supplementary Table 7) were also identified. No significant GO terms or pathways were identified for the 51 genes in the MP1-specific selective sweep regions. However, many genes were related to developmental process (Supplementary Table 8). In addition, several genes under MP1-specific selection were also related to cell or neuron projection as MP2 sub-population (Supplementary Table 8). Only SCN11A gene was commonly selected in MP1 and MP2 sub-populations (Supplementary Table 6). Selection signatures of this gene were not identified in other domestic pig populations. The SCN11A gene is involved in pain perception [44], suggesting that this gene has been specifically selected during minipig sub-populations to be adapted in laboratory environments.

A variance of genomic signatures related to body size in the minipig population

To determine whether there are common genetic mechanisms affecting body sizes of minipigs despite their large genetic variance, ZFst values were calculated for each minipig breed against the EDP population having the largest body size among domestic pigs (Methods). An average of 12,759 highly differentiated regions (ZFst > 2.5) containing an average of 1,582 genes were identified in each minipig breed (Supplementary Table 9). Significant GO terms and KEGG pathways for genes in selective sweep regions of each breed were very different, suggesting the diversity of selective sweeps in minipigs (Supplementary Fig. 8). The commonality of differentiated regions and genes was also divergent among minipigs (Fig. 4a). MP1 breeds shared relatively many differentiated regions, but only 305 genes were commonly differentiated among all MP1 breeds. Even in MP2 breeds, most of the differentiated regions and genes were breed-specific.

Fig. 4

Results of identifying differentiated genomic regions related to body size in the eight minipigs. a Left and right bar plots show the number of differentiated genomic regions and the number of genes in those regions in each minipig, respectively. Bar plots in middle panels represent the number of commonly differentiated genomic regions (top) and the number of common genes (bottom) shared among different combinations of minipigs. Colors of central circles indicate minipig sub-populations (green: MP1; yellow: MP2). b Candidate genes associated with body sizes of minipigs in GnRH secretion pathways. Candidate genes for different minipigs are denoted with different colors. The pathway was generated using pathview R package [45]. c,d Genotype plots for a specific 10 Kbp genomic region and its flanking region in the KISS1R (c) and GNA11 (d) genes highlighted with red boxes in b. Coordinates of the 10 Kbp genomic regions are shown as numbers on top lines. Yellow, dark blue, and red colors indicate homozygous reference, alternative genotype, and heterozygous genotype, respectively

Functional enrichment analysis for all genes located in differentiated genomic regions of minipig breeds altogether revealed many BP GO terms related to systemic, cell, and neural development (Supplementary Table 10). Many pathways in the KEGG analysis were associated with cell proliferation [46,47,48,49]. Selection patterns of genes involved in those GO terms and pathways were highly divergent among minipigs. For example, most genes in the GnRH secretion pathway were differentiated in at least one minipig compared to the EDP population (Fig. 4b). However, gene differentiation patterns were highly different among minipigs. In Bama minipig, cell surface receptors and ion channel genes including KISS1R (GPR54) were mainly differentiated from the EDP population (Fig. 4c). However, many signaling molecule genes such as GNA11 (Gq/11) were significantly differentiated in the Yucatan minipig (Fig. 4d).

Several other genes involved in body sizes of pigs, such as AR [8], LCORL [12], NR6A1 [12], and VRTN [13], were also observed in differentiated regions. For those genes, disparate differentiation patterns were also identified among minipigs (Supplementary Fig. 9). Compared to minipigs, allele patterns across EDP breeds were quite similar for all genes. These results suggest that biological mechanisms affecting body size are especially divergent in minipigs.

Characterization of candidate genes regulating body sizes of pigs

Although many genes related to body size were separately differentiated among minipigs, 71 genes were commonly differentiated against EDP breeds in all minipigs (Supplementary Table 11). To determine the direction of differentiation, $\pi$ ratio between the minipig and the EDP population was calculated (Methods). Genes in regions with the top and bottom 5% $\pi$ ratio values in both MP1 and MP2 sub-populations were considered as EDP- and MP-specific differentiated genes, respectively (Methods). Other genes were treated as genes differentiated in both populations. XP-EHH, Tajima’s D, and heterozygosity scores were also calculated for regions containing genes for MP1 and MP2 sub-populations against the EDP population (Methods).

Most genes were differentiated in both populations, whereas 26 genes and one gene were EDP- and MP-specifically differentiated, respectively (Fig. 5a, Supplementary Table 12). They included genes already known to be involved in body size, such as PLAG1 [12], CHM [50], and ESR1 [51], which showed distinct differentiation patterns in different populations (Fig. 5b). In the case of the ESR1 gene, which was found as an MP-specific gene, the gene contained two genomic regions with the bottom 5% $\pi$ ratio values in both MP1 and MP2 sub-populations (Fig. 5b, c). Especially, the region 14,465,001–14,475,000 on chromosome 1 showed much higher XP-EHH scores than genome-wide averages. This region also showed much lower heterozygosity and Tajima’s D scores for both MP1 and MP2 sub-populations than for the EDP population. Conversely, the overall PLAG1 gene body, which was deemed an EDP-specific gene, was in the genomic regions with the top 5% $\pi$ ratio values in both MP1 and MP2 sub-populations. These regions showed much lower XP-EHH scores than genome-wide average and lower heterozygosity and Tajima’s D scores in the EDP population compared to minipig sub-populations (Fig. 5d). The CHM showed differentiated patterns in both minipigs and EDP population. For example, although they showed high ZFst values, the $\pi$ ratio and XP-EHH values were close to the genome-wide average (Fig. 5e). Also, the heterozygosity and Tajima’s D scores were similar among populations.

Fig. 5

Statistics of candidate genes involved in body size that are commonly identified in all minipigs compared to European domestic pigs (EDP). a,b Nucleotide diversity (π) ratio distribution for MP1 and MP2 sub-populations against the EDP population for body size-related genes commonly identified in all minipigs. Each dot represents a single 10 Kbp genomic region within each gene. a Blue and purple dots indicate top and bottom 5% $\pi$ ratio values in both MP1 and MP2 populations. Other regions are represented by grey dots. b 10 Kbp genomic regions in PLAG1, CHM, and ESR1 genes are highlighted by yellow, green, and red color, respectively. c-e Distribution of Z-transformed fixation index (ZFst), XP-EHH, π ratio, Tajima’s D, and heterozygosity scores among different pig populations for ESR1 (c), PLAG1 (d) and CHM genes (e). Dotted horizontal lines with blue and orange color indicate genome-wide average of XP-EHH and π ratio scores for MP1 and MP2 sub-populations against the EDP population, respectively. Dotted vertical lines indicate the start and end position of each gene. The coordinate of region highlighted with grey color in (c) is 14.465–14.475 (Mbp) in chromosome 1

To interrogate the effect of selection on gene expression for the three genes, we conducted differential gene expression (DGE) analyses between minipigs and the EDP population in three tissues including brain, liver, and muscle. (Methods; Supplementary Table 13). The CHM gene was commonly down-regulated in the muscle tissue of minipigs compared to that in the EDP population, while the ESR1 gene was up-regulated in the brain tissue of minipigs compared to that in the EDP population. The PLAG1 gene did not show any consistent differential expression patterns between the two populations in any tissues.

Discussion

In this study, we performed a comprehensive population genetic analysis for eight minipig breeds using SNP data identified from a total of 41 pig breeds and five outgroup species covering various geographical regions and different body sizes. Using this large dataset, we revealed highly diverse and unique genetic structures within the minipig populations. They could be separated into two sub-populations, called MP1 and MP2 sub-populations, which were genetically close to Asian and European pigs, respectively. Admixture analyses showed that the MP1 inbred has a single Asian ancestry, whereas the MP1 outbred minipigs have both Asian and European ancestries, indicating admixture between the two ancestries. In contrast to MP1 minipigs, MP2 minipigs mainly have European ancestry with few Asian ancestry. We also found that the ancestry composition is highly variable even within the same minipig sub-population. Consistent with the result of admixture analyses, the MP1 sub-population shared a larger number of selective sweep regions during domestication with the Asian population than with the European population while the MP2 sub-population showed the opposite pattern.

Demographic signatures varied even within the minipig population. MP1 minipigs showed highly different ROH and π patterns between inbred and outbred minipigs. For example, the inbred Bama minipig showed the largest NROH and SROH and the lowest $\pi$ value, suggesting that Bama experienced a stronger population bottleneck than others. The Bama minipig also showed the largest frequency of long ROH ($\ge$ 8 Mbp), indicating that the population bottleneck has been maintained until recently. Although outbred MP1 minipigs showed smaller NROH and SROH values than inbred ones, their values were larger than those of non-admixed pigs. Thus, outbred MP1 minipigs also experienced weak population bottleneck, which might be caused by the founder effect during admixture events. MP2 minipigs showed highly similar $\pi$ values but highly different NROH and SROH values. Although the Minnesota minipig is an outbred one, NROH and SROH values were larger than those of the inbred Yucatan minipig. This might imply that the Minnesota minipig was developed by crossbreeding but was recently established as an inbred colony. The LD was more diverse among minipigs even in the same sub-population or those with similar breeding histories. These findings suggest that although minipigs share similar genetic ancestry or breeding history, they can have highly diverse demographic signatures.

We also identified several candidate genes involved in body sizes of pigs. As expected, these genes were related to systemic development and cell proliferation. Interestingly, those genes were also involved in neural development and the GnRH secretion pathway. The brain, an organ in the central nervous system, produces and secretes diverse growth hormones [52] and its size is correlated with body size [53]. Compared to other domestic pigs, many minipigs show early sexual maturation [54] which is regulated by GnRH [55]. Early sexual maturation has been reported to lead to short heights in human [56]. These findings suggest that body sizes of pigs are regulated by multiple mechanisms such as cell proliferation, systemic development, neural development, and hormone regulation. In addition, since several genes related to actin or cell projection important for cell growth and division [57] and neural development [58, 59] were identified in the minipig-specific selective sweep regions during domestication, the minipig-specific selection pressure affecting actin-related cell or neural development could exist.

Similar to demographic signatures, selection signatures within mechanisms regulating body size were also highly variable within minipigs. Regarding neural development, many brain development and neuron-related GO terms were specifically identified in the Minnesota minipig. However, several synapse-related GO terms were found in the Yucatan minipig. Some MP1 minipigs also showed different GO terms related to neural development. In addition, differentiation patterns of related genes in the GnRH secretion pathway were highly different among minipigs. We also found disparate selection signatures in several genes related to body size in pigs, such as NR6A1 [12], VRTN [13], AR [8], and LCORL [12]. Compared to the minipig population, the EDP population showed quite similar allele patterns with each other in those genes. This suggests that mechanisms regulating body sizes of minipigs differ by breed.

Despite the variability of genetic structures in minipigs, we found that 71 genes, including PLAG1, ESR1, and CHM, commonly differentiated in all minipigs against EDP breeds. PLAG1 is already known to be involved in pig body size [12]. ESR1 has been reported to be strongly associated with human height [51] and is related to the regulation of bone growth and maturation [60]. In the case of CHM, its disruption is known to inhibit normal development of the embryo and reduce body size in mouse [61] and zebrafish [62]. Therefore, these three genes can be candidate key genes regulating body sizes of different pig populations. In differentiation direction analysis, PLAG1 and ESR1 were differentiated only in a single population (EDP and minipig, respectively), while CHM was differentiated in both populations with different directions. When we compared the expression of those genes between minipigs and EDP breeds, CHM and ESR1 were highly differentially expressed in muscle and brain tissues, respectively, whereas PLAG1 showed similar expression levels in all tissues. These results suggest that allelic differentiation in CHM and ESR1 may cause changes in gene expression, whereas that of PLAG1 may be affected by different mechanisms such as differential splicing.

Conclusions

We found several unique and variable characteristics of genetic structures within the minipig population compared to other pig populations. Moreover, although minipigs have similar genetic structures, they show highly different demographic and selective signatures. Our results also suggest that those genetic variances could independently shape common phenotypes in the population, such as dwarfism, through highly different underlying genetic mechanisms. Our research findings could help elucidate biological mechanisms underlying development and provide a basis for establishing sustainable breed development program.