The data consisted of five commercial pig populations; Synthetic (Large-White based), Pietrain, Landrace, Large-White and Duroc ( Table 1), obtained from herds owned by Topigs Norsvin. The phenotypic data was composed by backfat records, measured at the end of the test period (average live weight of 120 kg) in each evaluated population. The data were classified in two datasets: PHENOTYPED, which consisted of all phenotyped animals and their contemporary’s animals; and GENOTYPED, which is a subset of phenotyped animals with genotypes and phenotypes (Table 1).

Using PHENOTYPED dataset, phenotypes were pre-corrected for all non-genetic effects to use as response variable in further analysis. For that, the non-genetic effects were estimated within population using the following animal linear model in ASReml v3.0 (Gilmour et al., 2009): BF ijkl = μ + sex i + hym j + a k + litter l + e ijkl (1) where BF ijkl is the phenotypic observations (BF) of the k animal; μ is the overall mean; sex i is the fixed effect of sex i; hym j is the fixed effect of the herd-year-month j of birth; a k is the random additive genetic effect of animal k; litter l is the random effect of litter l and e ijkl is the random residual effect. It was assumed that a ∼N (0, A σ a 2 ), l ∼N (0, I σ l 2 ) and e ∼N (0, I σ e 2 ), where σ a 2 , σ l 2 and σ e 2 are the additive genetic variance, litter variance and residual variance, respectively; A is a pedigree-based numerator relationship matrix and I is the identity matrix.

In total, 134,887 animals were genotyped across the different populations and SNP chip densities (Table 2). The GENOTYPED dataset includes animals genotyped using different Illumina’s medium density SNP chips (50, 60 and 80K), as well as animals genotyped with the high-density Axiom™ Porcine Genotyping Array 660K (Thermo Fisher Scientific, Waltham, Massachusetts, United States) (Table 2).

Genotype quality control analysis was performed within population and SNP chip to exclude SNPs with genotype call rate <0.95, minor allele frequency <0.01, strong deviation from Hardy-Weinberg equilibrium p-value < 1 × 10^-12, SNPs located on sex chromosomes and unmapped SNPs. The positions of the SNPs are based on the Sscrofa11.1 assembly of the reference genome (Warr et al., 2020a). All genotyped animals had a frequency of missing genotypes <0.05 and were therefore all kept for further analyses.

After the quality control, the medium density genotypes were imputed towards the high density genotypes within population using Fimpute v2.2 (Sargolzaei et al., 2014). After imputation, a total of 31,183; 16,189; 38,298; 36,414 and 11,274 animals and 488,585, 504,525, 443,619, 518,659 and 417,532 SNPs were available for the Synthetic, Pietrain, Landrace, Large-White and Duroc populations, respectively, and used for the GWAS.

A single-SNP GWAS was performed with the imputed GENOTYPED dataset for each population using the following linear animal model in GCTA software (Yang et al., 2011, 2014): y ∗ k = μ + X β ^ + u k + e k (2) where y ∗ k is the pre-corrected phenotype of the k animal (pre-corrected for all non-genetic effects); μ is the average of the pre-corrected phenotype; X is the genotype, coded as 0, 1, or 2 copies of one of the alleles of the k animal for the evaluated SNP; β ^ is the unknown allele substitution effect of the evaluated SNP; u_k is the residual polygenic effect, assuming u ∼N (0, G σ u 2 ), which accounted for the (co)variances between animals due to relationships by formation of an G matrix (genomic numerator relationship matrix build using the imputed genotypes), σ u 2 is the additive genetic variance; and e k is the random residual effect which was assumed to be distributed as ∼N (0, I σ e 2 ).

The proportion of phenotypic variance explained by a SNP was defined as σ S N P 2 σ P 2 , where σ P 2 is total phenotypic variance (sum of the additive and residual variances) which was estimated based on model (2) without a SNP effect, and σ S N P 2 = 2 p q ( β ) 2 , where p and q are the allele frequencies and β ^ the estimated allele substitution effect of the evaluated SNP. Variance components were also estimated using model (2) without a SNP effect and the heritability was defined as the proportion of σ P 2 explained by the genetic variance. The association between a SNP and the phenotype was declared significant when–log₁₀(p-value) > 8.

To search for possible causal mutations at the chromosome 5 backfat locus, fine mapping analyses using high-throughput sequencing (WGS, ChIP-Seq, ATAC-Seq and RNA-Seq data) was performed using following steps: (1) building and calculation of haplotype frequencies; (2) sequence data processing and evaluation of variants within the haplotype region; (3) evaluation of protein interactions and gene expression.

From the GWAS results, the most significant SNP (lead SNP) was identified. Subsequently we performed a haplotype analysis up-and downstream of the lead SNP. For that, firstly the Beagle v5.1 software (Browning et al., 2018) was used to phase the imputed genotypes taking 0.5 Mb downstream and 0.5 Mb upstream of the lead SNP (SSC5:66,103,958bp) in the QTL region SSC5:66 ± 0.5 Mb. The lead SNP and the 40 markers with closest vicinity of the lead SNP (20 markers on each side) were selected and used in the HaploPi pipeline (Oliveira, 2019) to build and calculate the haplotype frequencies. After performing the haplotype analyses, the core haplotype (SSC5: 66094630-66140348) was used to perform a haplotype-based association analysis using the number of copies of the core haplotype. For that, the same mixed-model described for the single-SNP GWAS (model 2) was used, in which the SNP genotypes were replaced by the genotype of the core haplotype (0, 1 or 2).

Whole-genome sequence (WGS) data was used to access the complete set of variants within the haplotype region (SSC5:66094630-66140348). For this purpose, 659 animals (Synthetic = 187; Pietrain = 40; Landrace = 227; Large White = 205) were available (Table 2). Genomic DNA from ear biopsies was extracted using BioSprint DNA Kit (Qiagen) and its concentration and quality was measured using NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific). The data was sequenced on Illumina Hiseq 2000 150 bp paired-end reads. The reads were aligned to the Sscrofa11.1 (Warr et al., 2020b) using BWA-MEM v0.7.15 (Li and Durbin, 2009) with an average mappability of 96.51% and a sample coverage ranging from 6.6–22.7X (∼10X average). SAMtools v1.9 dedup function was used to remove PCR duplicates (Li et al., 2009). Variant calling was performed with Freebayes v1.1.0 with following settings: min-base-quality 10, min-alternate-fraction 0.2, haplotype-length 0 and min-alternate-count 2 (Garrison and Marth, 2012). Variants with Phred quality score <20, and within 3 bp of an indel were discarded (Li et al., 2009). Variants were annotated using the Ensembl variant effect predictor (VEP, release 99) (McLaren et al., 2016). The impact of missense variants was predicted using the SIFT tool v6.2.1 (Kumar et al., 2009).

Chromatin immunoprecipitation–coupled sequencing (ChIP-Seq) marked by trimethylated histone H3 lysine 4 (H3K4me3), mono-methylation histone H3 lysine 4 (H3K4me1) and acetylated histone 3 lysine 27(H3K27ac) and ATAC-Seq data were downloaded from the European Nucleotide Archive (ENA) with accession numbers PRJEB28147 (ChIP-Seq - brain, liver, muscle and testis), and PRJNA665194 (ATAC-Seq – adipose, cerebellum, cerebral cortex, hypothalamus, liver, lung, skeletal muscle and spleen). ChIP-Seq marks and ATAC-Seq data were aligned against the pig reference genome (Sscrofa 11.1) using BWA-mem v. 0.7.17 (Li, 2013), pre-processed using SAMtools v. 1.13 (Li and Durbin, 2009) and broad and narrow peaks calling was performed on MACS v2.2.6 (Zhang et al., 2008) (parameters used for mapping, pre-processing and peak calling can be found in Supplementary Material S1.

Four different tissues: liver, spleen, lung and muscle were subjected to RNA extraction from 100 animals. The animals were F2 crosses resulting from a F1 sow (Landrace*LargeWhite) and a Synthetic boar line. The collected tissue was stored with RNAlater (Thermo Fisher Scientific) at -80°C until further use. RNA was extracted from the 400 samples using the QIAshredder homogenizer kit (Qiagen) and RNA extraction was performed with the Rneasy kit (Qiagen) following manufacture’s guidelines. Once all the samples passed the required quality control parameters, they were sequenced at 150 bp paired-end in an Illumina 6000 sequencing platform. DNA from 87 of these animals was also extracted from the spleen tissue and was subjected to high density genotyping with the Axiom™ Porcine Genotyping Array (Thermo Fisher Scientific) that queries 660K variant markers.

For the bioinformatic analysis, RNA-Seq reads were trimmed using Trim Galore v0.3.7 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). We mapped the RNA-seq data to the Sscrofa11.1 reference genome and Ensembl version 104 using STAR v2.7.8 (Dobin et al., 2013). Alignments with an alignment MAPQ score <30 were filtered using SAMtools v0.1.19 (Li and Durbin, 2009). Gene counts were determined using htseq-count v0.11.1, and then TMM-normalized as counts per million (CPM) with EdgeR (Robinson et al., 2009). Plotting of expression levels per genotype class was performed using the seaborn python package (Waskom, 2021). For the eQTL analysis, genotypes were firstly filtered as in (Derks et al., 2021) with plink v.1.9 (Purcell et al., 2007). Only genes ±1 Mbp from the target lead SNP of interest were used for the analysis. The single-SNP association analysis was carried for each of the tissues with the GCTA v1.25.3 software (Yang et al., 2011) with the following model: Y i = µ + S N P i + e i (3) where ( Y i ) is the CPM gene abundance modeled as a function of the population mean (µ), fixed effect of each SNP ( S N P i ), and a random residual effect ( e i ). The genetic variance explained by a SNP ( σ S N P 2 = 2pqα ² ) was estimated based on the allele frequencies (p and q) and the estimated allele substitution effect (α). The proportion of phenotypic variance explained by the SNP was defined as σ S N P 2 / σ P 2 , where σ P 2 is total phenotypic variance (sum of the additive and residual variances) which was estimated based on model (3) without a SNP effect. Significant association between a SNP and expression was detected using a p-value < 0.05.

In this study, the use of large-scale genotype data allowed us to identify a A/G SNP (rs80985094: 5:66,103,958) on SSC5 with a strong association with backfat across populations. This SNP showed significant association with backfat in four (Synthetic, Pietrain, Landrace and Large White) out of the five evaluated pig populations (Figure 1). The association between rs80985094 and backfat was not identified in the Duroc population because the G allele, associated with increased backfat is fixed.

The most significant SNP (rs80985094) showed a -log10 (p-value) of 57.67, 43.29, 67.61 and 37.91 in the Synthetic, Pietrain, Landrace and Large-White population, respectively. The lead SNP, located at 66,103,958 bp on SSC5, explains up to 6% of the genetic variance in backfat (Table 4). Among the four populations which showed significant association between the lead SNP and backfat, the frequency of the allele associated with increased backfat (G) is the highest in Large-White population (0.79). In the other populations, the G allele is the minor allele (Table 4).

We phased the genotypes using Beagle 5.1 (Browning and Browning, 2007) and selected the most frequent haplotypes either wich contains the A and G allele in the four populations. Based on the haplotype analyses, we identified a core haplotype with four SNPs surrounding the lead SNP (GCAAG) present in all populations (size ∼28 Kb). Assuming that recombination reduced the common interval surrounding a new variant reducing back fat, the region containing the causal mutation was expanded by two more SNPs (one SNP upstream and one SNP downstream) to search for possible functional variants. This increases the final region of interest to 46 kb to search for the causal variant (AGCAAGC = Synthetic; AGCAAGG = Pietrain; GGCAAGA = Landrace; AGCAAGG = Large-White). The frequencies of these core haplotypes are 57.70%, 35.87%, 22.67% and 16.85% for each population (Synthetic, Pietrain, Landrace and Large-White, respectively) (Table 5). No common core haplotype surrounding the lead SNP could be identified across the populations surrounding the G allele (Table 6). Hence, we believe that the G to A mutation happened on the GCGAG haplotype and was subsequently selected to produce pigs with less backfat. This underlines the hypothesis that the G allele is the wild-type variant and the A allele marks the derived new variant. In line with this hypothesis, the core haplotype containing the wild type G allele could be identified at low frequency. Linkage disequilibrium (r ²) between the variants in the core haplotype and the lead SNP (A) range from r ² = 0.03 and r ² = 0.96 (Table 7). The additive effect of the A and G allele on backfat across the evaluated pig populations is shown in Figure 2.

A haplotype-based association analysis was performed to compare the association with the lead SNP and the core haplotype (GCAAG) as defined in Table 5 (orange haplotype). The association with the core haplotype was lower than the lead SNP (Table 8). This gives further support that the lead SNP is causal.

For fine mapping, we evaluated a total of 659 whole genome sequenced animals, 222 sequenced animals are homozygous for the A allele, 252 animals are heterozygous, and 185 animals are homozygous for the G allele. Analyzing the data of these subset of animals, we found two variants (besides the lead SNP) with high frequency across populations and in high LD with the lead SNP (Table 9).

The LD between these variants and the lead SNP varies among the populations (Table 9). The lead SNP and SSC5:66097445 are annotated within the intronic region of the Cyclin D2 (CCND2) gene (ENSSSCG00000038694), while SSC5:66190273 is annotated in intergenic region upstream of the CCND2 gene.

We studied the expression of the CCND2 and neighboring genes in a population of 100 crossbred animals in four different tissues (liver, spleen, lung, muscle). Table 10 shows the association of the lead SNP with the expression of genes in the backfat locus (SSC5:65-67 Mb) in four tissues.

The results show that the expression of the CCND2 gene is significantly associated with the lead SNP in liver, lung and spleen (Figure 3). Also, the expression of the AKAP3 (all four tissues, although with low expression CPM: 3-5) and the C12orf4 (liver, muscle) gene show significant association with the lead SNP. The AKAP3 gene encodes a member of A-kinase anchoring proteins (AKAPs). This protein is highly expressed in spermatozoa may function as a regulator of motility, capacitation, and the acrosome reaction. The C12orf4 gene encodes a protein involved in mast cell degranulation.

The core haplotype region ranges from SSC5:66094630 to SSC5:66140348 and it covers the CCND2 gene. Therefore, to discover potentially regulatory elements, we used ChIP-seq and ATAC-Seq public datasets to investigate the region of interest. ChIP-Seq data was used to detect enhancers and promoters marked by H3K4me3, H3K4ac27 and H3K4me1. ATAC-Seq was used to detect open chromatin.

We have called broad peaks for marks on ATAC-Seq and ChIP-Seq (Figure 4 and Supplementary Material S2). All markers have peak signals just downstream of the lead SNP and are partially covered within the core haplotype region (Supplementary Text S1). In human a predicted promoter is located at the lead SNP location, further supporting a potential regulatory role of the lead SNP (Supplementary Material S3).