Genome-wide association and evolutionary analyses reveal the formation of swine facial wrinkles in Chinese Erhualian pigs

Descriptive statistics of phenotypic traits

In this study, we measured facial wrinkles using two different methods in a herd of Erhualian pigs involving 332 individuals. Descriptive statistics for both facial wrinkle phenotypes were shown in Table 1. The averages of wrinkle pixel ratio and wrinkle score were 0.09 ± 0.03 and 3.20 ± 1.04, respectively. The wrinkle pixel ratio ranged from 0.03 to 0.15, while the wrinkle score from 1 to 5. The two wrinkle phenotypes varied widely with the coefficients of 30.2% and 32.6%, respectively, which indicated that large divergence existed in the Erhualian pig population for the traits of facial wrinkles. And both phenotypes of facial wrinkles shared close relationship with a correlation coefficient of 0.777 (P < 0.01).

To estimate sex effect on facial wrinkles, we classed Erhualian pigs into male and female groups. The average of the wrinkle pixel ratio were 0.08 and 0.09 in males and females, while the average of the wrinkle score were 3.06 and 3.35 in males and females, respectively (Supplementary Figure 2). The folding extent of facial wrinkles in females was slightly heavier than males (P = 1.87×10^-3 and P = 9.93×10^-3). Thus the gender slightly affects facial wrinkles in the Erhualian population, and sows generally have more facial wrinkles than boars.

Heritability estimation

The heritability estimates of facial wrinkles were high, 0.68 for wrinkle pixel ratio and 0.74 for wrinkle score (Table 1), which indicated that the trait of facial wrinkles in the Erhualian pigs was highly heritable and might be determined by a remarkable genetic effect.

Single-trait GWAS and multi-trait meta-analysis

In the 60K SNP data, no individual was removed due to low genotyping rate of less than 0.05, while 955 SNPs were excluded for high missing rate of more than 10% and 24,897 SNPs were removed for low minor allele frequency of less than 0.01. Finally, 35,977 SNPs and 332 individuals were kept for the next association analyses.

Single-trait genome-wide association analyses were performed for both wrinkle phenotypes in the Erhualian pig population. To assess the plausibly of our observed data, we performed the “Q-Q” plot analysis for both wrinkle traits. The results showed that the distribution of observed P-values deviated from expected P-values in the end (Supplementary Figure 3). And the inflation factors (λ) of the GWAS for both traits are 1.017 and 1.023, suggesting that our GWAS results were reasonable. In the single-trait GWAS, we identified only one major QTL at 34.8 Mb on SSC7 for swine facial wrinkles, which contained 13 and 14 significant SNPs associated with wrinkle pixel ratio (Figure 1A) and wrinkle score (Figure 1B), respectively. The most significant SNP for both traits was same, rs80983858, which located at the 3,255 bp downstream of gene GMR4. It explained 21.80% of phenotypic variation for wrinkle pixel ratio and 31.3% for wrinkle score (Table 2).

The two traits used for facial wrinkles measurements are highly correlated (Table 1), but they had different focuses: the pixel ratio mainly focused at the heavy wrinkle points, and the wrinkle score laid emphasis on the overall appearance. To validate the results of the single-trait GWASs and possibly increase the QTL power, we performed a two-trait meta-analysis. The two-trait meta-analysis was conducted based on the SNP effects from the single-trait GWAS. The QTL at 34.8 Mb on SSC7 was convinced with 13 significant SNPs passing over the suggestive significance threshold (1.39× 10⁻⁶), in which the most significant SNP was also rs80983858 with P = 4.45 × 10^-13 (Figure 1C and Table 2). The effect of G allele is benefit to increase facial wrinkles (Figure 1D).

Genome-wide gene association analysis

Genome-wide gene association analyses were performed for two phenotypes of facial wrinkles using MAGMA software in the Erhualian pig population. We identified only one major significant associated region at 34.8 Mb on SSC7 (Supplementary Figure 4), which was consistent with the result of SNPs-based GWAS. The top 10 significant genes for facial wrinkle score were GRM4, DST, ZNF451, MLIP, CDKAL1, PRIM2, PRL, RNF144B, LOC100513868 and FAM83B (Supplementary Figure 4A); and the top 10 genes for wrinkle pixel ratio were ZNF451, GRM4, MLIP, PRIM2, REC114, PRL, DST, TTC6, CDKAL1 and STXBP6 (Supplementary Figure 4B). Among these genes, GRM4, DST, ZNF451 and MLIP were strong candidate genes affecting the formation of swine facial wrinkles.

Evolutionary and selection analyses

To determine whether the candidate region affecting swine facial wrinkles was under selection or not, we have carried out three analyses on this region, including Tajima’s D test, genetic differentiation estimate of F_ST and calculation of XPEHH score. First, we merged the present 60K SNP dataset of 332 Erhualian pigs with previous 60K SNP dataset of 21 Chinese wild boars [9, 10]. We performed Tajima’s D test on the candidate region for Erhualian pigs and Chinese wild boars. Low Tajima’s D values less than -1 were detected in the Erhualian pigs but not in Chinese wild boars (Figure 2A), which suggested that this region was under selection in Erhualian pigs. We also found the region centered on 35.8 Mb ranging from 35.1 Mb to 36.8 Mb was highly differentiated between Erhualian pigs and Wild boars (Figure 2B). And high XPEHH scores were detected in this candidate region (Figure 2C). To validate the above results, we further performed these analyses on this region using our previous whole-genome sequencing dataset of Chinese Erhualian pigs [11, 12] and Chinese wild boars [12]. On the candidate region from 34.7 Mb to 36.2 Mb, the Tajima’s D values less than -1 were detected in the Erhualian pigs but not in Chinese wild boars (Supplementary Figure 5A). Meanwhile, high F_ST values (more than 0.15) and large XPEHH scores (more than 0.9) between Erhualian pigs and Chinese wild boars were observed on the target candidate region (Supplementary Figure 5B and 5C). These results strongly supported that the candidate region affecting the formation of facial wrinkles was under positive selection.

In addition, we found that the haplotypes containing G allele, benefiting to wrinkle formation in Erhualian pigs, were clustered into a single branch (Figure 2D). This result further suggested that the haplotypes with favorable wrinkle formation were under selection in Chinese Erhualian pigs.

Ethics statement

All procedures used for this study and involving animals are in compliance with guidelines for the care and utility of experimental animals established by the Ministry of Agriculture of China. The ethics committee of Jiangxi Agricultural University specifically approved this study.

Experimental animals

A herd of Chinese Erhualian pigs (n=332) were used as experimental animals as described previously [33]. Briefly, the Erhualian pigs at the age of 80 ± 3 days from nine sires and 52 dams were bought from the national conservation farm for the Erhualian breed in Changzhou, Jiangsu province, which covered all sire lines and most of dam lines raised in the conservation farm. All piglets were raised twice per day with a corn- and soybean-based diet containing 16% crude protein, 3100 kJ digestive energy and 0.78% lysine at a farm in Nanchang, Jiangxi province, and with free access to water. Male pigs were castrated at the age of 60 days. All pigs including 168 males and 164 females were uniformly slaughtered at 300 ± 3 days at an abattoir for trait measurements.

Phenotype measurements

After slaughter, swine heads were cut off at the joint between the first cervical vertebra and occipital bone. Front head picture of each individual was taken under same light condition using the camera of Xiaomi mobile phone Mi 3. We applied two methods to estimate the folding extent of facial wrinkles for all the experimental animals.

For the first method, wrinkle pixel ratio was calculated by comparing the sketch area of facial wrinkles to the whole area of facial contour with Photoshop software (Supplementary Figure 6A). First, each photo of the front head was imported into Photoshop software and was trimmed to remain the whole face by crop tool. Then we regulated the size of trimmed figure to fit a blank image with width size of 4 cm (472 pixels) and height size of 6 cm (709 pixels) without changing the width-to-length ratio by rectangular marquee tool and free transform tool. The images of 64 individuals arranged by 8 × 8 were placed on a 32 × 48 cm² canvas, and then were printed on A4 paper with default option of "Scale to Fit Media". Their facial wrinkles and peripheral facial contours were traced on translucent copy paper over the printed images using HB pencil. The traced wrinkles and peripheral facial contours were scanned and reimported into Photoshop software. The area pixels of the traced facial wrinkles (TFW) and blank facial contours (BFC) were calculated by Photoshop software. Wrinkle pixel ratio was considered as TFW/BFC, representing one phenotype for the folding extent of facial wrinkles.

For the secondary method, we established a scoring criteria for swine facial wrinkles from 1 to 5, in which 1 represents minimum folding extent, 5 represents maximum folding extent (Supplementary Figure 6B). Three volunteers were convened to score the facial wrinkles in all experimental pigs at the scoring step of 0.5 according to the scoring criteria. The three volunteers were trained to read the score according to the scoring criteria of swine facial wrinkles before their formal scoring. If the score difference for a same sample was larger than 1 among the three independent volunteers, this sample would be selected out and independently rescored based on the scoring criteria by three volunteers until the difference was equal to or smaller than 1. Finally, the average values of wrinkle scores determined by these three volunteers were calculated and treated as another phenotype for the folding extent of facial wrinkles.

SNP genotyping and quality control

Genomic DNA was extracted from ear tissues using a standard phenol/chloroform method. DNA concentration and quality were measured by a Nanodrop100 spectrophotometer (Thermo Fisher, USA) and agarose gel electrophoresis. All DNA samples were diluted to a final concentration of 50 ng/μl with ultraviolet absorption ratios of 260/280 > 1.8 and 260/230 > 2.0. Then these samples were genotyped for the Porcine SNP60 Beadchips on an iScan System (Illumina, USA) following the manufacturer’s protocol. Quality control of SNP genotype data was performed by PLINK (version 1.07). SNPs with a call rate < 90% and minor allele frequency (MAF) < 1% were removed from further analyses; individuals with missing genotypes > 10% or Mendelian errors > 5% were also excluded. Finally, 35,977 SNPs, including 477 SNPs from the X-chromosome were kept for GWAS analyses. Since male individual contains only one X chromosome, SNPs on the X-chromosome for males were treated as homozygous genotypes for further association analyses.

Heritability estimates

Heritabilities of the measured traits were estimated based on the genotype data using the polygenic function implemented in R package of GenABEL [34]. A maximize polygenic model [35] was used as follows:

$\text{logL}=\left(\text{1}/\text{2}\right)\left[{\text{log}}_{\text{e}}\left|\text{Σ}\right|+{\left(\text{y}-\text{μ}\right)}^{\text{T}}{\text{Σ}}^{-\text{1}}\left(\text{y}-\text{μ}\right)\right],$

in which y indicates standardized phenotypic value, μ stands for the average standardized phenotypic value, Σ denotes the variance-covariance matrix, and Σ = AV_A + IV_R. Here, A is the relationship matrix of kinship, which is estimated from genetic data; V_A is the additive genetic variance due to polygenes; I stands for the identity matrix; and V_R is the residual variance. Heritability estimates were calculated with the equation h² = V_A/(V_A + V_R).

Single-trait genome-wide association studies

Association between each SNP and phenotypes of facial wrinkles in the population was evaluated under a generalized linear mixed model using polygenic and mmscore function implemented in GenABEL package [34]. The statistic model is shown as follows:

$y=X\beta +S\alpha +Z\mu +e,$

in which y is the vector of phenotypes; β is the estimator of fixed effects including sex and batch; α is the vector of SNP substitution effects; and μ is the random additive genetic effect following a multivariate normal distribution μ ~ N(0, GV_G), in which G is the comparability kinship matrix between each individual estimated by whole-genome SNP as described in [36] and V_G is the polygenic additive variance; X, S, and Z are the incidence vectors (matrixes) for β, α, and μ, respectively; and e is a vector of residual effects normally distributed N(0, lV_E), V_E is the residual variance.

Genome-wide significance thresholds were set as P_0.05 = 0.05/N with Bonferroni adjustment, and the suggestive significance threshold was set to be 1/N, where N is the number of filtered SNPs in the data set. In this study, the thresholds were set as 1.39 × 10^-6 (0.05/35977) and 2.78 × 10^-5 (1/35977), respectively. The 2-LOD drop-off method was used to determine confidence intervals of the identified loci [37]. Percentage of variance (%) explained by each top significant SNP in the GWAS was calculated using the following formula:

$\text{var}=\left[\left({\text{MS}}_{\text{reduce1}}-{\text{MS}}_{\text{full}}\right)/{\text{MS}}_{\text{reduce}}\right]\times \text{1}00,$

in which MS_reduce1, MS_full, and MS_reduce indicate the mean square (MS) in the linear models including 3 effects (mean, sex and batch), 4 effects (mean, sex, batch and SNP), and only mean, respectively.

Multi-trait meta-analysis

A Chi-square statistic was calculated for the importance of the effects of SNP_i (i=1, 2, 3, …, 35977) across both wrinkle phenotypes estimated in the present study, which approximately follows a Chi square distribution with 2 degrees of freedom. For each SNP, the multi-trait statistic was calculated by the formula [38]:

${\chi }_{multi-trait}^2=t_i^{\text{'}}{V}^{-1}{t}_i$

where t_i is a 2×1 vector of the signed t values of the ith SNP for the two wrinkle traits, t’_i is the transpose of vector t_i(1×2), V^-1 is the inverse of 2×2 correlation matrix between two traits is the correlation over the 35,977 SNPs effects(signed t values) of the two traits.

Genome-wide gene association analysis

The SNP-based P values from the single SNP-by-SNP GWAS were used as input for the gene-based analysis. We used all 16,852 Sus scrofa protein-coding genes from the NCBI gene definitions as the basis for a genome-wide gene association analysis (GWGAS) in MAGMA (http://ctg.cncr.nl/software/magma). After SNP annotation, there were 5,520 genes that were covered by at least one SNP. Gene association tests were performed taking LD between SNPs into account. We applied a stringent Bonferroni correction to account for multiple testing, setting the genome-wide threshold for significance at 9.06 × 10⁻⁶.

Evolutionary analyses for the target region

The target region centering around the QTL from 32 to 38 Mb on SSC7 was chose for the evolutionary and selection analyses employing two datasets. The SNP data of 21 Chinese wild boars [9, 10] were included and merged together with the present data of Erhualian pigs to form the first 60K SNP dataset. Previous whole- genome sequencing data of 19 Erhualian pigs [11, 12] and 6 Chinese wild boars [12] with average depth of ~25X were combined together to form the second sequencing dataset. Genetic differentiation estimates of F_ST between Erhualian and wild boar populations were calculated as described in Weir’s paper [39] using the SNP dataset on the target region. Tajima’s D within both Erhualian and wild boar populations were estimated as described in Tajima’s paper [40]. The average F_ST and Tajima’s D were calculated in sliding windows of 200 kb with step size of 100 kb. Selection footprints on the target region were also detected by a between-population method of XPEHH implemented in Selscan software [41].

Phases of the 25 SNPs within the target region (from 34,634,619 to 35,422,882 bp on SSC7) were reconstructed using Phase v2.1.1 [42] with default options. Phased alleles were linked together to form haplotype sequences, which were used for reconstruction of neighbor-joining (NJ) phylogenetic trees by MEGA7 [43] with 1,000 times bootstrapping.