Guangdong Provincial Key Lab of Agro-Animal Genomics and Molecular Breeding/National Engineering Research Center for Breeding Swine Industry/College of Animal Science,South China Agricultural University,Guangzhou 510642,P.R.China

Abstract In modern pig breeding programs,growth and fatness are vital economic traits that significantly influence porcine production.To identify underlying variants and candidate genes associated with growth and fatness traits,a total of 1 067 genotyped Duroc pigs with de-regressed estimated breeding values (DEBV) records were analyzed in a genome wide association study (GWAS) by using a single marker regression model. In total,28 potential single nucleotide polymorphisms (SNPs)were associated with these traits of interest. Moreover,VPS4B,PHLPP1,and some other genes were highlighted as functionally plausible candidate genes that compose the underlying genetic architecture of porcine growth and fatness traits. Our findings contribute to a better understanding of the genetic architectures underlying swine growth and fatness traits that can be potentially used in pig breeding programs.

Keywords:pig,GWAS,growth trait,fatness trait,candidate gene

1.Introduction

Growth and fatness traits are vital economic traits that significantly influence porcine production in pig breeding programs. Porcine production is mainly measured by days to 100 kg (D100),loin muscle area and backfat thickness.Estimated heritability of the D100 trait ranges from 0.31 to 0.47 (Li and Kennedy 1994;Tanget al.2019),and there are 103 quantitative trait locus (QTLs) for the D100 trait in the Pig Quantitative Trait Locus Database (QTLdb) in the Animal QTLdb (https://www.animalgenome.org/cgi-bin/QTLdb/index,release 37). The loin muscle area also has a moderate-heritability between 0.45 and 0.50 (Suzukiet al.2005;Nakanoet al.2015),and 327 QTLs have been identified being associated with this trait. There are 450 QTLs for the backfat thickness trait in the Pig QTLdb,and these QTLs are distributed across the entire porcine genome(Qiaoet al.2015;Guoet al.2017b;Jianget al.2018).

With the advance of the high-density single nucleotide polymorphism (SNP) array,genome-wide association study(GWAS) has been widely applied to detect the variants associated with traits of interest in high resolution. GWAS has been used to explore the underlying genetic architecture of various traits,such as reproduction (Chenet al.2019),growth and fatness traits (Guoet al.2017a) in pigs.Additionally,the results of GWAS can play a guiding role in manual selection and can cover the shortage of traditional artificial selection.

Although a number of genomic regions and genes associated with growth and fatness traits have been identified in previous studies,more genes underlying these complex traits are yet to be uncovered. Furthermore,the previously reported genes need to be confirmed in different populationsviaGWAS. We hypothesized that some previously reported and novel genes associated with porcine growth and fatness traits would be identified in an independent population. Hence,we conducted a GWAS in a Duroc population to identify QTL regions associated with two growth traits,D100 and loin muscle area at 100 kg (L100),and one fatness trait,backfat thickness at 100 kg (B100).Then,functional enrichment to the identified QTL regions and comparison between the results of current study and previous studies were carried out. This study offers useful information for pig breeding programs.

2.Material and methods

2.1.Population and phenotyping

The research population was normally maintained in a breeding herd in south China (Fujian,China). Two growth traits,D100 and L100,and one fatness trait,B100,were calculated with default parameters and recorded for all pigs using the herdsman swine management software(Appendix A).

Records on the three traits of 4 539 Duroc pigs from this population were collected from 2008 to 2017 for further investigation. A multi-traits animal model was used to estimate heritability and genetic correlation. The model is as follows:

where Y is the vector of phenotypic records,b is the vector of fixed effects,including sex and year-season,a is the vector of additive genetics,pe is the vector of permanent environmental effects,e is a vector of residuals,and X,Z,and W are incidence matrices for b,a,and pe,respectively.

The estimated breeding values (EBVs) of all pigs and the reliabilities of EBVs were imputed using animal model best linear unbiased prediction (Henderson 1975). Then,de-regressed EBV (DEBV) were calculated for each pig to remove the contribution of information from relatives(VanRaden and Wiggans 1991). The equation is as follows:

where DEBV is the de-regressed EBV,PA represents the parental average,and EBV and REL are the estimated breeding value and reliability of each pig,respectively.

2.2.Genomic DNA extraction and genotyping

Genomic DNA was extracted from the ear tissues of 1 067 Duroc pigs (81 boars and 986 sows) using the TaKaRa MiniBEST Universal Genomic DNA Extraction Kit (Ver 4.0).The OD260/280ratios of DNA samples were quantified with NanoDrop 200 (Thermo Scientific,USA). The samples with an OD260/280ratio between 1.7 and 2.0 were genotyped using either the Illumina PorcineSNP60 BeadChip (Illumina,San Diego,CA,USA) containing 61 565 SNPs or the GeneSeek GGP-Porcine chip (Neogen Corporation,Lansing,MI,USA)containing 50 697 SNPs. There are 33 359 common SNPs in these two chips,and the missing genotypes were imputed using BEAGLE Software (Browning and Browning 2007).

2.3.Data quality control

Quality control for genotypes was conducted with the criterion below:firstly,individuals with a call rate ＜90%were excluded. Then,all SNPs were filtered according to the following criteria:(1) call rate＜90%;(2) minor allele frequency (MAF)＜1%;(3) Hardy-Weinberg Equilibrium(HWE) testingP-value＜10–6. This step was performed using PLINK Software (Purcellet al.2007). After quality control,32 147 SNPs were remained for further analysis.The descriptive statistical analysis for phenotypes was performed using R language (R Foundation for Statistical Computing 2018).

2.4.Statistical analysis

The GWAS was performed with the following single marker regression model by GEMMA Software (Zhou and Stephens 2012):

where y is a vector of dependent variable (DEBVs in this study),μ is the overall mean,G is the realized relationship matrix constructed with markers,I is an identity matrix,is the additive genetic variance,Z represents incidence matrices corresponding to u,e is the vector of residual errors,andis the residual variance.

To confirm the thresholds for the genome-wide significance and suggestive significance,the number of effectively independent tests were calculated and reported in a previous study (Zhanget al.2019). The genomewide significance threshold used in the current study was 0.05/9266=5.40×10–6,and the genome-wide suggestive significance threshold was 1/9266=1.08×10–4.

2.5.Functional enrichment analysis

Candidate genes were identified according to their physical positions and functions based on theSus scrofa10.2 reference genome assembly (http://may2017.archive.ensembl.org/index.html). SNP-containing annotated genes for each potential SNP were obtained and identified as candidate genes. Genes that located on less than 1 Mb away from the potential SNPs were selected for functional enrichment analysis. Gene ontology (GO) terms annotation was conducted using the R package clusterProfiler (Yuet al.2012),and the threshold by the Benjamini-Hochberg method (adjustedP-value＜0.05). The extended regions of the potential SNPs that overlapped with the reported QTL regions were considered to be signals associated with the signals of a given trait,while the reported QTLs were retrieved from the Pig QTLdb within the Animal QTLdb (Huet al.2019).

3.Results

3.1.Description of phenotypes

Descriptive statistics of DEBVs for the traits of D100,L100,and B100 are presented in Table 1. In the present study,the heritability of D100,L100,and B100 was 0.300,0.380,and 0.215,respectively.

3.2.Genome-wide association results

Manhattan plots and Q–Q plots across the whole genome of the three growth and fatness traits are shown in Fig.1 and Appendix B,respectively. The inflation factors λ of D100,L100,and B100 was 1.02,0.99,and 0.96,respectively,which were deemed to be safe. Additionally,the genomewide significant SNPs and the SNPs reaching the suggestive significance are listed in Table 2.

For the D100 trait,a total of 12 potential SNPs was identified,andVPS4B,PHLPP1,ENSSSCG00000034988,CDH20,ENSSSCG00000004911,andENSSSCG00000033637were identified as SNP-containing annotated genes.WU_10.2_2_80933825,WU_10.2_2_83165289,and WU_10.2_6_150933 were found to be associated with the L100 trait,andADAMTS2was shown to be an SNPcontaining gene. Furthermore,13 SNPs were identified being potentially associated with the B100 trait,andNUDT3,RPS10,TSHZ1,PHLPP1,ENSSSCG00000034988,CDH20,ENSSSCG00000004911,TXN2,GRM4,andSNRPCwere found to be SNP-containing annotated genes. Interestingly,there were five common SNPs and four common SNP-containing annotated genes,ranging from 17.63 to 17.93 Mb,on chromosome 1 for the D100 and B100 traits.

3.3.Functional enrichment results

To annotate the potential SNPs,the candidate genes overlapping with the extended genomic regions were selected for GO terms enrichment analysis (Fig.2). For the D100 trait,115 candidate genes were identified. These genes were enriched in certain GO terms related to enzyme regulator or inhibitor activity. Additionally,a total of 184 candidate genes were found to be associated with the B100 trait,while these genes were shown to be involved with metabolism-associated GO terms. Although 111 candidate genes overlapped with the extended genomic regions of the potential SNPs associated with L100,no GO terms met the threshold.

Then,the information of Pig QTLdb was used to further annotate the extended genomic regions. The top 10 traits with the highest enrichments are shown in Fig.3. For the D100 trait,the average daily gain QTL showed the highest enrichments,with loin muscle area and backfat traits being enriched as well. For the L100 trait,the loin muscle area QTL was among the top 10 QTLs with the highest enrichments. Additionally,various backfat-associated QTLs were highly enriched for the B100 trait. The average daily gain and loin muscle area QTL were also highly enriched.

4.Discussion

We performed a GWAS in a Duroc population to explore the genetic architectures of growth and fatness traits. A total of 28 significant or suggestive SNPs and several candidate genes were identified for the D100,L100,and B100 traits. These results hence confirmed our hypothesis that some important variants and genes potentially associated with porcine growth and fatness traits have not yet to be identified. In this study,some novel genes have been identified to be significantly associated with porcine growth and fatness traits,suggesting that previously reported geneswere far from complete. Therefore,it is necessary to perform GWAS in different pig populations to identify more genes underlying the complex traits of growth and fatness,which is beneficial for pig breeding programs.

Table 1 Phenotypes (de-regressed estimated breeding values (DEBV)) of three porcine growth and fatness traits1)

Fig.1 Manhattan plots of three porcine growth and fatness traits distinguished by text labels. The y-axis of Manhattan plots displays the–log10 (P-value) of each single-nucleotide polymorphism (SNP) in the genome-wide association analysis. The two horizontal lines divided SNP with P-value＜5.40×10–6 and ＜1.08×10–4,respectively. The red dots stand for the potential SNPs associated with the traits of days to 100 kg (D100),loin muscle area at 100 kg (L100),and backfat thickness at 100 kg (B100).

Previous studies concluded that some SNP-containing annotated genes were highly associated with the three porcine growth and fatness traits. For the D100 trait,six SNP-containing annotated genes were identified.Vacuolar protein sorting 4 homolog B (VPS4B) involves in the endosomal sorting complexes,and is necessary for transport machinery (Tanget al.2015). VPS4B has been reported to be crucial for the degradation of membrane receptors,regulation of epidermal growth factor receptor and insulin receptor (Liuet al.2013). For the L100 trait,ADAM metallopeptidase with thrombospondin type 1 motif 2 (ADAMTS2) is the only identified SNP-containing annotated gene.ADAMTS2,a protein coding gene,was found being associated with severe skin fragility and joint laxity (Porteret al.2005). Ten SNP-containing annotated genes were found for the B100 trait. Among them,nudix hydrolase 3 (NUDT3),ribosomal protein S10 (RPS10) and glutamate metabotropic receptor 4 (GRM4) have been reported being related to porcine growth and carcass traits in previous studies (Wanget al.2014;Gonget al.2019;Zhuanget al.2019). Furthermore,GRM4has been identified as a candidate gene that increases facial wrinkles in Chinese Erhualian pigs (Huanget al.2019) and affects feed efficiency in purebred castrated Yorkshire pigs (Houet al.2018). Teashirt zinc finger homeobox 1 (TSHZ1),another candidate gene,is likely to influence the marbling rate evaluation trait in Nellore cattle (Fonsecaet al.2020).Small nuclear ribonucleoprotein polypeptide C (SNRPC)has been identified as a signal in exploring the genetic architecture of human thinness compared to severe obesity(Riveros-McKayet al.2019).

Interestingly,it was found that the D100 and B100 traits share four common SNP-containing annotated genes. PH domain leucine-rich repeat-containing protein phosphatase 1(PHLPP1) encodes a phosphatase,and has been reported as a functionally plausible candidate gene associated with average daily weight,backfat thickness,body weight,and carcass weight in pigs (Howardet al.2015;Guoet al.2017b). Cadherin 20 (CDH20),involved in cell adhesion pathway (Linet al.2014),was reported being associated with porcine growth and fatness traits for the first time. We also identifiedENSSSCG00000034988,a long non-coding RNA (lncRNA),as a candidate gene associated with these traits for the first time. Since there were some studies showing that variants in lncRNA could similarly influence phenotype of complex traits (Zhenget al.2016;Jinet al.2017;Zhanget al.2018),ENSSSCG00000034988deserves further investigations.ENSSSCG00000004911is a novel O-acyltransferase like protein-like gene,and has been identified as a candidate gene associated with the traits of daily feed intake and daily feeder occupation time in pigs(Reyeret al.2017).

Table 2 Potential single-nucleotide polymorphisms (SNPs) and candidate genes detected in the genome-wide association study for three porcine growth and fatness traits

Fig.2 The enrichment results of days to 100 kg (D100) and backfat thickness at 100 kg (B100) traits distinguished by text labels.

Fig.3 The enrichment results of three porcine growth and fatness traits,days to 100 kg (D100),loin muscle area at 100 kg (L100),backfat thickness at 100 kg (B100),based on the Pig Quantitative Trait Loci (QTL) Database (https://www.animalgenome.org/cgi-bin/QTLdb/SS/index). Each dot denotes a list of the extended regions overlapped QTLs. The top 10 traits with the highest enrichments are shown using red triangles.

The results of the potential SNPs and candidate genes were confirmed through the enrichments of the three porcine growth and fatness traits based on the Pig QTLdb and GO terms. These findings help to better understand the genetic architectures of these three traits and benefit pig breeding programs.

5.Conclusion

Through a genome-wide association study on three growth and fatness traits in a Duroc pig population,we detected 28 potential SNPs that were associated with these traits of interest.VPS4B,PHLPP1,and some other genes were highlighted as functionally plausible candidate genes that compose the underlying genetic architecture of porcine growth and fatness traits. These findings could provide more information on the genetic architectures of porcine growth and fatness traits and have the potential to be used in pig breeding programs.

Acknowledgements

This research was supported by the earmarked fund for China Agriculture Research System (CARS-35),the National Natural Science Foundation of China (31772556),and the Key R&D Program of Guangdong Province,China(2018B020203002).

Declaration of competing interest

The authors declare that they have no conflict of interest.

Ethical approval

Animal care and experiments were conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology,China,revised in June,2004) and were approved by the Animal Care and Use Committee of South China Agricultural University,Guangzhou,China (permit number:SCAU#2013-10).

Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm