Whole-genome re-sequencing association study for direct genetic effects and social genetic effects of six growth traits in Large White pigs

Socially affected traits are affected by direct genetic effects (DGE) and social genetic effects (SGE). DGE and SGE of an individual directly quantify the genetic influence of its own phenotypes and the phenotypes of other individuals, respectively. In the current study, a total of 3,276 Large White pigs from different pens were used, and each pen contained 10 piglets. DGE and SGE were estimated for six socially affected traits, and then a GWAS was conducted to identify SNPs associated with DGE and SGE. Based on the whole-genome re-sequencing, 40 Large White pigs were genotyped and 10,501,384 high quality SNPs were retained for single-locus and multi-locus GWAS. For single-locus GWAS, a total of 54 SNPs associated with DGE and 33 SNPs with SGE exceeded the threshold (P < 5.00E-07) were detected for six growth traits. Of these, 22 SNPs with pleiotropic effects were shared by DGE and SGE. For multi-locus GWAS, a total of 72 and 110 putative QTNs were detected for DGE and SGE, respectively. Of these, 5 SNPs with pleiotropic effects were shared by DGE and SGE. It is noteworthy that 2 SNPs (SSC8: 16438396 for DGE and SSC17: 9697454 for SGE) were detected in single-locus and multi-locus GWAS. Furthermore, 15 positional candidate genes shared by SGE and DGE were identified because of their roles in behaviour, health and disease. Identification of genetic variants and candidate genes for DGE and SGE for socially affected traits will provide a new insight to understand the genetic architecture of socially affected traits in pigs.

the pen effects and showed that the heritability of SGE was approximately were about 0.003 11 and 0.0007 12 for ADG. The reported social heritabilities were 0.007, 0.001, and 0.003 for body weight, carcass backfat depth, and carcass lean meat content in heavy pigs 13 , respectively. Moreover, Rostellato et al. (2015) reported that the social genetic effect contributed largely to total genetic variance for the carcass traits (ranging from 33.2% to 35%) in pigs 13 . Canario et al. (2017) showed that the heritability of SGE was 0.001 to 0.004 for pig growth rate and its genetic correlations among direct and SGE were −0.28 to 0.12 in a Swedish pig population 3 .
In general, the pigs were grouped into different pens, and the social effects among members strongly affected pig productive performance. In pigs, the estimation of the social genetic effect was mainly focused on ADG 14 , feed intake, backfat and muscle depth 7 . Several models that considered different fixed factors were built to estimate SGE. The animal model was first extended to SGE in Japanese quail 15 . This model incorporated the competitive effects in classical BLUP (best linear unbiased prediction), and substantially resulted in an extra selection response. Muir 16 and Bijma et al. 17,18 presented a full model that simultaneously contained both direct genetic effects and social genetic effects. They estimated SGEs for growth traits and found that SGE contributed the main part of genetic variance in growth rate and feed intake in pigs. In summary, the previous studies on SGE mainly focused on the development of an evaluation model and the estimation of genetic parameters. The genetic architecture and causal pathway for SGE are still unknown.
Pig growth traits are affected by direct genetic effects (DGE) and SGE. To investigate the genetic architecture of DGE and SGE, this study conducted a GWAS for DGE and SGE of six growth traits using the whole-genome re-sequenced data from the Large White pig population. The purpose of this study was to estimate the DGE and SGE for six socially affected traits, and reveal the genomic variation and candidate genes for DGE and SGE of six growth traits in pigs.

Results
Summary of estimated DGE and SGE. Using this full model Y = Xb + Z D a D + Z S a S + Wl + Vg + e (see Materials and Methods), the DGE and SGE were estimated for each socially affected trait in Large White pigs. Table S1 listed the describe statistics about the DGE and SGE for these six traits. The DGE and SGE contained residual effects of 40 genotyped pigs were separately used to conduct GWAS. Variance components were presented in Table 1 for six socially affected traits.
The single-locus GWAS results using GEMMA. For DGE. A single marker test was performed to identify SNPs associated with DGE of 6 growth traits. The Table 2 and Fig. 1 show the GWAS results for DGE of 6 growth traits. And Q-Q plots are shown in Supplementary Fig. 1. A total of 54 SNPs and 20 genes were detected for DGE in this study.
For SGE. To explore the potential genetic variants associated with SGE, a whole-genome re-sequencing association study was conducted for SGE of 6 growth traits. All of the Manhattan plots are shown in Fig. 1 and the Q-Q plots are shown in Supplementary Fig. 2. At the genome-wide threshold (P < 5 × 10 −8 ), 4 identified SNPs www.nature.com/scientificreports www.nature.com/scientificreports/ were associated with SGE of days to 100 kg (sgeD100), 7 with sgeRFI. At the lenient threshold (P < 5 × 10 −7 ) for suggestive associations, 1 SNP was associated with SGE of average daily gain (sgeADG), 14 with sgeD100, 2 with sgeADFI, 4 with sgeFCR, and 1 with sgeRFI, but no association was found for sgeB100 (Table 3 and Fig. 2).
A total of 8 SNPs were found on SSC1 and SSC16 for sgeRFI (Table 3 and Fig. 2F). Two SNPs located on SSC1 and five SNPs located on SSC16 were found to reach the genome-wide significance threshold (P < 5 × 10 −8 ). The top SNPs (SSC16: 43805766) with the lowest P-value of 3.72 × 10 −9 showed the strongest association effect in the  www.nature.com/scientificreports www.nature.com/scientificreports/ A total of 15 significant SNPs were shared by DGE and SGE that were distributed on SSC1, SSC2, SSC3, SSC4, SSC7, SSC9, SSC14 and SSC16 (Table 4). In the region of SSC14: 98.69-98.73 Mb, the top SNP (SSC14: 98711408 bp) was located in the PRKG1 gene for dgeADG and sgeADG. In the region of SSC2: 14.29-14.33 Mb, the top SNP (SSC2: 14306048 bp) was located in the LOC110259270 gene, and in the region of SSC8:112. 41-112.45 Mb, the top SNP (SSC8: 112431517 bp) was located in GAR1 and the CFI gene for dgeD100 and sgeD100, respectively.
The multi-locus GWAS results using FarmCPU. The summary of GWAS results using the FarmCPU approach is listed in Tables S2 and S3. In total, the current study identified 51 and 101 SNPs for DGE and SGE, respectively. For DGE, 4 SNPs were associated with dgeB100, 45 with dgeD100 and 2 with dgeFCR (Table S2), respectively. For SGE, 3 SNPs were associated with sgeADG, 63 with sgeB100, 18 with sgeD100, 7 with sgeADFI, 5 with sgeFCR and 4 with sgeRFI (Table S3), respectively. Among these SNPs, a total of 7 SNPs were shared by dgeD100 and sgeD100 (Table 4). Additionally, three candidate genes were found at 20 kb around these top SNPs for dgeD100 and sgeD100. Among theme, two genes, NR3C2 and PKD2, were separately located in SSC8: 80.56-80.60 Mb and 131.00-131.04 Mb. One gene, LYPLAL1, was located in the region of SSC10: 8.96-9.00 Mb.
The multi-locus GWAS results using FASTmrEMMA. The (Table S5), respectively. In the region of SSC17: 9.68-9.72 Mb, a QTN located on SSC17: 9697454 bp was identified by single-locus GWAS and multi-locus GWAS for sgeD100. This QTN with 6.04 QTN effect explained 14.30% of the phenotypic variance. Furthermore, a total of 5 putative QTNs were shared by SGE and DGE of each trait (Table 4). For the ADFI, B100 and D100 trait, the QTN SSC10: 11015344 was shared by DGE and SGE. Four QTNs including SSC4: 140729546 bp, SSC8: 33116699 bp, SSC11: 15329359 bp and SSC13: 21796761 bp, were detected in both dgeD100 and sgeD100. Manhattan plots of genome-wide association analysis results associated to SGE for six growth traits using GEMMA software. The threshold for genome-wide and suggestive significance was set at P = 5.00 × 10 −8 and P = 5.00 × 10 −7 , respectively.  www.nature.com/scientificreports www.nature.com/scientificreports/ Interestingly common loci in the three methods. This study performed a GWAS for DGE and SGE using the GEMMA, FarmCPU and FASTmrEMMA methods. Interestingly, two identified SNPs were validated in three methods (Table 5). One SNP (SSC8: 16438396 bp) was associated with dgeB100. One SNP (SSC17: 9697454 bp) was associated with sgeD100.

Discussion
The estimation of genetic effect. Many studies have reported that the social genetic effect are widespread in hens 19 , mice 20 , and pigs 7 . However, at present, very few GWAS for SGE are conducted in pigs because the genetic parameters of SGE are difficult to measure. In past several decades, the heritability of social effects were estimated to be very low in pigs 11,21,22 using the common REML method. Therefore, social effects were often ignored in past pig breeding. Instead of common heritability of social effects, the ratio between total heritable variance and phenotypic variance was used as a measure of the importance of SGE. On that basis, the social effect was found to contribute the heritable variance in growth rate and feed intake 23 . This study used the full model (including DGE and SGE) described by Bergsma et al. 23 to estimate DGE and SGE of six growth traits. In this model, on one hand, the fixed effects contained sex, tested year and month, the random litter effect and the random group effect were corrected. On the other hand, additional 3236 phenotypic records for ADG, D100, B100, ADFI, FCR and RFI and pedigree information were added into mixed model equations, which improved the prediction accuracy of DGE and SGE for these traits.
The genetic architecture for SGE. Currently, the genomic information has increased the understanding of complex traits, but the majority of studies have focused only on the DGE of studied traits. Published studies demonstrated that ignoring SGE may severely bias estimates of DGE, and the genetic basis of SGE may be different from the genetic basis of DGE 8 . If SGE contributes to the studied traits, the study of genetic architecture should simultaneously consider DGE and SGE for socially-affected traits. However, few studies have focused on the genetic architecture of SGE. This study simultaneously quantified the SGE and DGE, and performed GWAS to reveal the    24 identified several genomic regions and candidate genes associated with SGE for ADG in Landrace pigs. In their study, the individual ADG and the average ADG of unrelated pen-mates were directly used as DGE and SGE, respectively. In our study, the SGE and DGE were estimated for six growth traits, using the full animal model. Then, (1) pseudo phenotype values including DGE and residual effects were used to perform GWAS for DGE of six growth traits; and (2) for SGE of six growth traits, pseudo phenotype values including SGE and residual effects were used to perform a GWAS. Using these pseudo phenotype values, this study revealed the genetic architecture of DGE and SGE for six growth traits in Large White pigs. This study detected 151 and 205 SNPs for DGE and SGE, respectively. These results did not confirm the results with Hong et al. 24 . In the present study, the power of QTL detection could be limited by the low estimated heritability of SGE and the small population size in this study. Furthermore the assessment of SGE was different between our study and the previous study 24 . In a previous study, the average ADG of unrelated pen mates was used to detect genetic variants of SGE. However, in our study, the full model Y = Xb + Z D a D + Z S a S + Wl + Vg + e (see Materials and Methods) were used to estimate the SGE of each trait, and then the estimated SGE wasused to perform GWAS for SGE of each trait. Our study offered a new insight and effective method for the genetic architecture of socially affected traits in pigs.

Comparing the single-locus and multi-locus GWAS.
To minimize the number of false positives and obtain credible results, the SNPs with MAF > 0.1 were used to perform single-locus and multi-locus GWAS using GEMMA 25 , FarmCPU 26 and FASTmrEMMA 27,28 software. Due to the multiple tests, the single-locus GWAS required Bonferroni correction 29 . However, the quantitative traits are controlled by numerous polygenes with large or minor effects, and the Bonferroni correction is overly conservative to identify important loci and may result in false negative results 27 for complex traits. Recently, some multi-locus GWAS were proposed to detect significant SNPs for complex traits 30,31 . The significant threshold of the multi-locus GWAS is less stringent than single-locus GWAS. Because of the multi-locus and shrinkage natures of this method, a less stringent selection criterion is used in the multi-locus GWAS 31 . Although the multi-locus GWAS used a less stringent criterion, the study used simulated and real data to demonstrated that the multi-locus GWAS using the FASTmrEMMA model is more powerful and has less bias of QTN effect estimation than other methods 31,32 . The FASTmrEMMA method with high power, low false-positive rate and low computing time in QTN detection were used to perform multi-locus GWAS. Using GEMMA, 54 and 33 SNPs reached the threshold of 5.00 × 10 −7 for DGE and SGE, respectively. Using FarmCPU, 51 and 101 SNPs reached the threshold of 5.00 × 10 −7 for DGE and SGE, respectively. For FASTmrEMMA, a total of 46 and 71 putative QTNs were detected for DGE and SGE, respectively. Clearly, different methods provided different results. However, all these results show that numerous identified SNPs are associated with pleiotropic effects for multiple socially affected traits. Importantly, two SNPs were detected in three methods. Among them, one SNP (SSC8: 16438396 bp) was located in the region 16.42-16.46 Mb for dgeB100. In this region, four reported QTLs were associated with the B100 trait in pigs 33,34 . One SNP (SSC17: 9697454 bp) was detected in the region of 9.68-9.72 Mb for sgeD100. A previous study reported that numerous QTLs were associated with body weight 35 , behaviour 36 , immune capacity 37 and disease 38 in pigs. This SNP may contribute to weight, behaviour, immune and disease traits, and further affect socially affected traits. Thus this study suggested that the SNP SSC17: 9697454 bp was an informative causal locus for SGE in socially affected traits.
The promising candidate genes identified by GWAS. To investigate the genomic features that contribute to SGE and DGE in pigs, 40 Large White pigs were re-sequenced and 6 socially affected traits were studied. Based on the whole-genome re-sequencing data, this GWAS revealed the genetic architecture for 6 socially affected traits. The previous study only reported a total of 5 QTLs related to socially affected traits in pigs 39 . To date, no genes with their biological functions were reported to associate with SGE in animals. Thus, it is essential to reveal the candidate genes and understand the genetic architecture of DGE and SGE for future pig breeding programmes. In the current study, using a GWAS study, many candidate genes associated with DGE and SGE were found. A total of 20 genes for DGE and 9 genes for SGE were found using single-locus GWAS. A total of 42 and 54 genes were found for DGE and SGE using multi-locus GWAS, respectively. All of these genes were first reported to be associated with SGE in pigs. Additionally, a total of 15 candidate genes associated with pleiotropism for socially affected traits. Among them, 4 genes for single-locus GWAS (including PRKG1, uncharacterized LOC110259270, GABRR2 and ADAMTS6 genes) and 11 genes for multi-locus GWAS (including SLC26A2, HMGXB3, SPTB, PLEKHG3, SLC30A9, NRSN2, FOXO1, LYPLAL1, NR3C2, PKD2 and TRIB3 genes) were shared by DGE and SGE, which may imply the pleiotropism of these shared genes for socially affected traits in pigs. These candidate genes were play a role in behaviour, disease and health.
The PRKG1 gene is a protein coding gene that acts as a key mediator of the cyclic guanosine monophosphate (cGMP) signalling pathway and plays an important role in the cellular signal transduction process. In platelets, hippocampal neurons, smooth muscle and cerebellar Purkinje cells, the PRKG1 gene was found to be strongly expressed 40,41 . This gene is known to be associated with behaviour 41 , circadian rhythms 42 , sleep 42 , learning 42 and memory 41 . In humans, the PRKG1 plays an important role in stress response related traits 43 .
SORCS3 is an orphan receptor of the VPS10P domain receptor family, a group of sorting and signalling receptors central to many pathways in the control of neuronal viability and function. SORCS3 plays an important role in the nervous system. This gene was associated with behavioural activities in mice 44 .
On the region of SSC1: 57.25-57.29 Mb, GABRR2 was found to be associated with dgeFCR and sgeFCR in pigs. GABRR2 is a protein coding gene that encodes a receptor for gamma-aminobutyric acid (GABA). The GABA receptor regulates the neurotransmitter in the brain 45 , and plays an important role in the behavioral stress response and physiology in animals and humans [45][46][47][48] . In mice, the GABA was associated with the aggressiveness and sociability towards conspecifics 47 . The performance and physical condition were affected by GABA in heat-stressed Roman hens 49 . (2019) 9:9667 | https://doi.org/10.1038/s41598-019-45919-0 www.nature.com/scientificreports www.nature.com/scientificreports/ The strongest association (SSC16: 43805766, P = 4.85 × 10 −9 ) was located in a disintegrin and metalloproteinase with thrombospondin motifs 6 (ADAMTS6). The ADAMTS6 gene is a member of the ADAMTS protein family and is regulated by the cytokine TNF-alpha 50 . The functions of this gene were still not clearly elucidated and mainly reported about aetiology and played a role in the turnover of the extracellular matrix 51 . ADAMTS6 was expressed differently in each tissue, but the accurate function and the substrates of ADAMTS6 protein were not clearly demonstrated by a previous study. In past decades, numerous association studies have found the ADAMTS6 gene to be associated with complex traits in humans. Inguinal hernias 52 and central corneal thickness and keratoconus 53 were both identified ADAMTS6 in a GWAS, which suggested that the ADAMTS6 gene would affected the collagen homoeostasis in tissues and disorders and lung function. A previous study also suggested that this gene was related to osteosarcoma 54 . Some studies also demonstrated that the ADAMTS6 gene was significantly upregulated in cancer cells and stromal cells 55 . Importantly, this gene affected the intelligence and growth development, and the balanced translocation disruption of the ADAMTS6 gene was resulted in short stature and intellectual disability 56 . Given previous studies, ADAMTS6 gene was associated with growth and disease in biological organisms. The ADAMTS6 gene was first identified to associate with the sgeRFI. The sgeRFI was not only affected the feed efficiency, but also affected the growth and health of livestock. Considering the pleiotropic effect of this gene, the ADAMTS6 gene with a potential effect on socially affected traits would be identified as a prominent candidate gene.
In summary, numerous SNPs and genes were identified for DGE and SGE of six traits using single-locus and multi-locus GWAS in Large White pigs. These findings were particularly interesting to better understand the genetic and physiologic mechanisms of both DGE and SGE. Although, a limited number of individuals were used in this study, this study provided a new insight to investigate socially affected traits. Further study using a large population size should contribute to validating the genetic mechanisms for socially affected traits in pigs.

Materials and Methods
Ethics statement. All experimental procedures were performed in accordance with the Institutional Review Board (IRB14044) and the Institutional Animal Care and Use Committee of the Sichuan Agricultural University under permit number DKY-B20140302.
Animals and Trait Recorded. The phenotypic data was collected in 2017, including ADG, days to 100 kg (D100), backfat thickness to 100 kg (B100), average daily feed intake (ADFI), residual feed intake (RFI) and feed conversion ratio (FCR). Initially, the tested pigs were selected between 120 and 130 day of age, with the body weight (BW) about 60 kg. Then a total of 40 Large White pigs from 4 pens were used for this research, containing 27 female and 13 male pigs, and each pen contained 10 piglets homogenous in body weight and age in this study. In this test, these pigs were grouped in the standard commercial pens, which were fed ad libitum by an automatic feeder (The feed intake recording equipment of OSBOREN). The test began at about 65 kg (BW 1 ) and the ended at about 110 kg (BW 2 ) in the feeding trial.
During the test period, we measured the phenotype data every ten days, and calculated the average value in finally. The ADG, ADFI and FCR were directly calculated from the collected data. The D100 and B100 were calculated as below (Kennedy and Chesnais, unpublished data): The average metabolic body weight (AMW) was calculated for each individual using the following classical formula 58 : . .

AMW
To investigate the direct and social genetic effects, the following full model was built for each trait in this study.

D D S S
where Y is the vector of phenotypic observations; b is the vector of fixed effects, including sex, tested year and month; a D is a vector of DGE; a S is a vector of SGE; l is vector of random litter effects; g is vector of random group effects; e is random residual vector; X, Z D , Z S , W and V are incidence matrix of b, a D , a S , l and g, respectively. The variance-covariance matrix of a D and a S is denoted as: where A is additive genetic correlation matrix; σ A 2 d , σ A 2 ds and σ A 2 s are genetic variance and covariance between direct effects and social effects. At the i th row of Z S , the group members of individual i were set to 1 and the others were set to 0. For improving the estimation accuracy of DGE and SGE, other 3,236 phenotypic records about ADG, D100, B100, ADFI, FCR and RFI that were obtained from previous performance test also were added into this analysis. These traits were analyzed using the average information restricted maximum likelihood (AI-REML) algorithm by DMU software 59 . Supplementary Table 1 lists the describe statistics about the DGE and SGE for these six traits.
Genotyping. The pig's ear tissue was collected from Large White pigs and stored in 75% alcohol. The genomic DNA was extracted by the standard phenol/chloroform method. The Nanodrop-2000 spectrophotometer was used to control the quality of genomic DNA. To obtain nucleotide polymorphism information, a total of 40 Large White pigs were re-sequenced using the Illumina HiSeq PE150 platform, then about 2,000 Gb sequence data was obtained in total. The average sequencing depth of these samples were close to 20×. All the sequence reads were filtered for data quality and mapped to the Sscrofa11.1 reference sequence using BWA software 60 . The mapped reads were realigned by GATK software 61 . A total of 21,104,245 variants were identified using GATK software 61 . In order to minimize the number of false positives in this limited population size, the SNP with MAF > 0.1 were retained for further analysis. Thus, the nucleotide variants were filter based on the quality requirement with minor allele frequency (MAF > 0.1), Missing rate (Miss < 0.1), Hardy-Weinberg equilibrium (HWE < 1.0 × 10 −6 ), and read depth (dp > 6) using VCFtools 4.2 62 . Then, the SNPs on the sex chromosome and scaffolds were removed. After the quality control, a total of 10,501,384 SNPs across 18 autosomes were analyzed for association study.
Single-locus association analysis. The genome-wide efficient mixed-model analysis (GEMMA) 25 with a univariate liner mixed model were used to perform the single-locus association for DGE and SGE. The analysis models was written as follows: where y is a vector of DGE or SGE contained residual effects from previous model (1), X is the incidence matrix of SNP effects, m is the vector of SNP effects, W is the incidence matrix of residual polygene effects, a is the vector of residual polygene effects, e is the random residual vector. Notably, where y is a vector that phenotype value was corrected by fixed effect, DGE (for SGE), SGE (for DGE), random litter and random group effects.
Recently, several statistical methods have been contributed to determine the significance threshold, such as Bonferroni correction 63 , False discovery rate 64 , and Sidak correction method 65 . In published literatures, the traditional Bonferroni correction method was routinely adopted, but this threshold was not absolute. When identified SNPs exceeded the accepted genome-wide statistical significance threshold P < 1.0 × 10 −8 , the GWAS results are most reliable 66,67 . The identified associations (5.00 × 10 −8 < P ≤ 5.00 × 10 −7 ) also could be replicated from subsequent studies 28 . The proper significant threshold is various due to different populations, different traits, and different genotypic data [68][69][70] . Furthermore, the Bonferroni correction is too conservative due to the ignorance of linkages between SNPs, especially for the sequencing data where many adjacent SNPs are highly linked. In the sequencing data, the number of markers used in the Bonferroni correction should be more relevant to the number of segments of chromosome. Thus, the whole-genome significance threshold 5.00 × 10 −8 was established 71,72 and the suggestive significant threshold was set at 5.00 × 10 −7 in this study 73 . Multi-locus association analysis. The multi-locus association analysis were performed by FarmCPU 26 and the fast multi-locus random-SNP-effect EMMA (FASTmrEMMA) 31,74 . FarmCPU iteratively used fixed effect model and random effect model. Using the FASTmrEMMA, this study performed a GWAS for DGE and SGE for