Comparison of SLAF-seq with genotyping BeadChip
In the present study, we performed a GWAS based on SLAF-seq to screen and select candidate SNPs for BFT, CLP and CFP traits on the four-way crossbred population with Landrace, Yorkshire, Duroc and Saba pigs as the hybrid parents. SLAF-seq as a high-resolution method for genome-wide genotyping is considered as an enhanced reduced representation sequencing, which has several significant advantages: accurate genotyping by depth paired-end sequencing, reduced sequencing costs, pre-design strate for fragmentation efficiency, and double barcode system for large populations [24]. In the study, a total of 1,552,377 SLAFs from the pig reference genome (Sscrofa 11.1_102) were predicted using RsaI and HaeIII enzyme combination. With an average sequencing depth of 11.94-fold, a total of 10,784,484 SNPs were obtained from 223 crossbreds using SLAF-seq and a Q30 value reached 90.74%. After filtering, 227, 921 SNPs distributed among 38 pig chromosomes genome were used to perform GWAS for BFT, CLP and CFP traits, which were nearly three times the number of SNPs on Illumina PorcineSNP80 Genotyping BeadChip. Furthermore, compared with Genotyping BeadChip, SLAF-seq can provide new variation information and is more suitable for the identification of new variation sites. Currently, SLAF-seq was successfully applied to pig genotyping and identified a large number of new mutation sites [31, 32]. Considering the huge differences in genomic variation between Chinese and Western pig breeds, our results showed that SLAF-seq was a powerful method and had great potential for further study in more pig breeds. Therefore, SLAF-seq method could be considered a more competitive choice in pig genome research.
Comparison of QTLs identified in this study with findings of previous studies
Currently, 70, 273, 58, 586 and 15 QTLs associated with SBFT, LRBFT, LBFT, ABFT and 67RBFT, respectively, deposited in the pig QTL database (http://www.animalgenome.org/ cgi-bin/QTLdb/index, Dec 27, 2021). In the study, a total of 14 SNPs were identified as significant for the BFT traits investigated. Among the detected SNPs, four, three, two, seven and five SNPs, were significantly associated with SBFT, LRBFT, LBFT, ABFT and 67RBFT, respectively. For BFT traits, SNPs identified in the study were mainly distributed in SSC3, SSC7, SSC8, SSC9, SSC12, SSC13 and SSC15, which showed that the genetic structure affecting BFT traits was complex. Among significant SNPs, one SNP (SSC7:30292654) associated with ABFT was located in the QTL region that has been reported in the White Duroc × Erhualian F2 population [7]. Besides, four shared SNPs related to at least two BFT traits, including SSC7:21392136, SSC7:21466553, SSC7:29486003 and SSC7:29503670 were also identified on SSC7. The four loci were located in the QTL region that has been reported in the five domesticated modern populations (Large White, Duroc-Large White synthetic, Yorkshire/Large White, and Landrace) [33]. In the study, a total of five SNPs significantly associated with BFTs were found on SSC7 (between 21.39Mb and 30.29Mb). A lot of previous studies have shown that QTLs related to BFT traits were located in the region near the five loci on SSC7 [8, 34–36].
Furthermore, the locus of SSC15:117022248 related to SBFT was located between microsatellites SW120 and SW906, which has a high density of QTLs associated with fat content [37]. The SNP at SSC3:34308396 associated with LRBFT and 67RBFT was located in the QTL region reported in the Duroc-Pietrain population [38], and SSC8:36710706 related to LBFT was located in the QTL region and has been reported in White Duroc × Erhualian intercross population [5]. SSC13:197751407 associated with ABFT was located in the QTL region (195.0-206.7Mb), and has been reported in the F2 Duroc × Pietrain population [39]. Two SNPs at SSC9:8763434 and SSC9:8763503 related to 67RBFT were located in the QTL area, which has been reported in Landrace [40] and Meishan × Large White populations [41]. The SNP at 50896936 bp on SSC12 associated with 67RBFT was located in the QTL region (23.7-60.0Mb) reported in the Duroc-Pietrain population [42].
For CLP and CFP traits, 265 and 58 QTLs have been accumulated respectively in the pig QTL database (http://www.animalgenome.org/cgi-bin/QTLdb/index, Dec 27, 2021). In the study, a total of seven SNPs were identified as significant for the CLP and CFP traits. Among the detected SNPs, six and four SNPs, were significantly related to CLP and CFP traits, respectively. For CLP and CFP traits, SNPs were mainly distributed in SSC9, SSC10, SSC11 and SSC14. The SNP at SSC10:46524739 associated with CLP and CFP traits located in the QTL region (46.08-46.98Mb) of feed conversion ratio that has been reported in Duroc pigs [43]. Besides, two SNPs at SSC11:47266057 and SSC11:47266119 associated with CLP and CFP traits located in the QTL region and SSC9:8744721 associated with CLP traits located in the QTL region were reported in F2 Duroc × Pietrain populations [44]. SSC14:38533013 related to CLP was located in the QTL region of average backfat thickness that has been reported in Chinese Meishan × Dutch Large White populations [45] and populations of backcrossing F1 Meishan-White composite females to either Meishan or White composite boars [35]. Finally, SSC10:50888840 associated with CFP was located in the QTL region of shoulder subcutaneous fat thickness and has been reported in Meishan × Large White populations [41].
Of the 21 significant SNPs associated with BFT, CLP and CFP traits found in four-way crossbred pigs in our study, except for two significant SNPs (unmapped to the porcine reference chromosome) that have not been reported in previous studies, other SNPs were located in previously reported QTLs. All significant SNPs can be used as new candidate loci for future studies of porcine BFT, CLP and CFP traits. However, more validations need to be conducted in more pig breeds.
Possible pleiotropic QTL
Our results indicated that some loci might have pleiotropic effects on different traits. Interestingly, A 74,417 Kb interval (from 21.39 to 21.46 Mb) on SSC7 was shown to contain two SNPs that were associated with three BFT traits, including SSC7:21392136 for SBFT (P = 2.31E-06) and ABFT (P = 1.01E-05), and SSC7:21466553 for LRBFT (P = 1.11E-06) and ABFT (P = 1.05E-06). There was also an interval with a length of 0.8Mb (from 29.48 to 30.29 Mb) on SSC7, which was shown to contain three SNPs that were associated with four related traits, including SSC7:29503670 for SBFT (P = 3.94E-06), LRBFT (P = 8.98E-06), LBFT (P = 1.02E-06) and ABFT (P = 3.26E-08), SSC7:29486003 for SBFT (P = 5.29E-06) and ABFT (P = 2.12E-07), and SSC7:30292654 for ABFT (P = 6.58E-06). Coincidently, Qiao et al. also reported several significant SNPs for BFT traits on SSC7 in the 2.1Mb region between 33.30 Mb and 35.42 Mb [7]. On SSC9, two SNPs including SSC9:8763434 and SSC9:8763503 were associated with 67RBFT (P = 8.00E-06, P = 9.62E-06), and were close to SNP SSC9:8744721 which was associated with CLP (P = 9.18E-06). The result of the correlation analysis showed that the correlation coefficient between 67RBFT and CLP was negative and highly significant (r = − 0.47, P < 0.001) (Fig. S2). Three adjacent SNPs associated with both traits were within a region of 0.19Mb (between 87.44Mb and 87.63Mb) on SSC9, which indicated that this region might contain a pleiotropic QTL (Table 2).
Furthermore, the correlation coefficient between CLP and CFP was negative and highly significant (r = − 0.86, P < 0.001) (Fig. S2). Interestingly, two contiguous SNPs on SSC11, including SSC11:47,266,057 and SSC11:47,266,119 were significantly associated with both traits. Besides, a 4.36 Mb region (between 46.52 and 50.88 Mb) on SSC10 was associated with CLP and CFP, which indicated that the region might contain a pleiotropic QTL (Table 2).
Candidate genes for BFT traits
The GWAS result showed that six SNPs were associated with at least two BFTs at the significance level (p < 1.0E-05). The adjacent and within genes of the SSC7:21392136 and SSC7:21466553 were ZNF184 and ZNF391. GO analysis showed that ZNF184 and ZNF391 participated in regulation of transcription, DNA-templated and nucleic acid binding, etc. As is known, zinc-finger proteins (ZFPs) represent the largest transcription factor family in mammals [46]. These two genes belong to ZFP transcription factors, which play a crucial role in regulating diverse growth and development processes through nucleic acid binding and transcription activation, etc. [47]. An increasing number of ZFPs involved in adipogenesis have been discovered, such as zinc finger protein 423 (Zfp423) [48, 49], Zfp467[50, 51], Zfp521[52], ZNF395[53]. A large number of ZFPs play roles in preadipocyte differentiation, such as ZNF638[54], GATA2[55], GATA3[55] and SLUG [56], which belong to the C2C2-type zinc finger protein subfamily consisting of a highly conserved zinc finger DNA binding domain. It is known that ZNF184 and ZNF391 also are the C2H2-type zinc finger transcriptional factors. According to their structure and function, it is inferred that ZNF184 and ZNF391 might be involved in regulating adipogenesis as C2H2-type transcriptional factors. In addition, SSC7:29486003 and SSC7:29503670 were located in the intron region of COL21A1. The results of the functional analysis showed that COL21A1 was enriched in collagen trimer and extracellular matrix (ECM). It has been reported that COL21A1 is VAdomain-containing collagen with a domain structure and which is a part and the smallest of the FACIT family of collagen expressed in various tissues including heart, stomach, placenta, skeletal muscle, kidney and liver [57, 58]. The co-expression of collagen XXI and collagen I in tissues and muscles has an important role in the organization of interstitial collagen fibrils by connecting them to other matrix components or cells [57, 59]. A study finds that adipocyte differentiation is regulated by the deposition of collagen which is the major component of ECM [60]. So, COL21A1 plays role in ECM architecture and function. The ECM remodeling is associated with the modulation of adipogenesis during adipose tissue expansion [61]. Besides, Tanaka et al. study the structure of pig genome region containing fat-related QTLs by constructing bacterial artificial chromosome (BAC) overlap from DST to SRPK1 region on SSC7. The analysis of nucleotide sequences of three BAC clones including rearrangement points shows that COL21A1 and DST are located near the rearrangement point, and they do not exist in the corresponding human region [62]. Our findings would help to further understand the role of COL21A1 in backfat metabolism and this gene should be considered a strong candidate gene for BFT traits.
Furthermore, four genes, including OR4P4L, OR4C13L, OR5D14, and OR5D13L were located within 100 kb upstream and downstream of AEMK02000449.1:179574, which was significantly associated with ABFT and 67RBFT. Interestingly, OR1J4L was found nearby the AEMK02000598.1:813208, which was significantly related to CLP. In the study, GO analysis showed that OR4P4L, OR4C13L and OR1J4L were mainly enriched in G-protein coupled receptor signaling pathway and G-protein coupled receptor activity, etc. These five genes belong to the olfactory receptor (OR) gene family. The ORs are G protein-coupled receptors that recognize a large number of organic compounds in accordance with their homologous ligand, which have been found to be expressed in different organs or tissues in mammals [63, 64]. The ORs expressed in extra-nasal tissues are involved in a variety of biological processes, including the regulation of glucose and lipid metabolism [65]. Therefore, it could be inferred that these five OR genes might be involved in regulating the energy metabolism and lipid metabolism of pigs to affect BFT and CLP traits of pigs. However, there were no genes found within ± 100kb surround of SSC3:34308396, which was significant associations with LRBFT and 67RBFT. As above, these seven significant SNPs and eight genes might be important for the fatness of the pigs and be potential targets of BFT traits for further molecular breeding research.
In the remaining significant SNPs, BARD1 and ENSSSCG00000048197, located near SSC15:117022248 were found. Enrichment analysis showed that the BARD1 gene was involved in RNA binding and ubiquitin-protein transferase activity. At present, these two genes were not been reported in fat-related studies. Four genes, including GRM4, SMIM29, NUDT3 and HMGA1, were located within 100 kb upstream and downstream of SSC7:30292654. Of these four genes, HMGA1 can be a candidate gene for the locus, as HMGA1 is functionally related to fat metabolism and its variants have been associated with backfat thickness [7]. Some studies show that HMGA1 can serve as a regulator of IGF1 activity to regulate glucose uptake [66] and can bind with PPARG, a key modulator of adipocyte differentiation and glucose homeostasis [67]. Besides, the adjacent and within genes of the SSC13:197751407 were MRPS6 and SLC5A3. MRPS6 encodes a subunit of mitochondrial ribosomal protein 28s23 [68], while SLC5A3 is a gene embedded in MRPS6, encoding a protein that transports sodium and inositol under hyperosmotic stress [69]. Currently, these two genes were not reported in fat metabolism for further study. Five genes, including LIPT2, KCNE3, PGM2L1, P4HA3 and ENSSSCG00000042110 were nearby the two adjacent SNPs (SSC9:8763434 and SSC9:8763503), which were significantly associated with 67RBFT. Interestingly, the five genes were also found near SSC9:8744721, which was significantly related to CLP. Among the five genes, PGM2L1 was the closest gene to the three loci. GO analysis showed that PGM2L1 was enriched in carbohydrate metabolic process, phosphorylation and glucose-1,6-bisphosphate synthase activity, etc. It is known that PGM2L1 belonging to a distinct family of α-phosphohexomutases widely distributed in prokaryotes can make glucose-1,6-bisphosphate, a key metabolic regulator [70]. PGM2L1 plays an indispensable role in glycogen metabolism, glycolysis and gluconeogenesis as a key enzyme [70]. Our findings provided a new understanding of the role of PGM2L1 in the backfat metabolism of pigs and should be considered a strong candidate gene for BFT and CLP traits. Finally, the adjacent and within genes of the SSC12:50896936 were PITPNM3, PIMREG and AIPL1. PITPNM3 (phosphatidylinositol transfer protein 3) and AIPL1 (Aryl hydrocarbon receptor-interacting protein-like 1) are previously reported in association with vision (i.e., phototransduction [71, 72]). The PIMREG gene may be related to porcine immunity [73]. Alvarenga et al. find that the PITPNM3 gene located on BTA19 and SSC12 is associated with temperament in cattle and struggling bouts (frequency and duration) and latency at an age of 12 days in pigs, which is a key gene related to animal behavioral characteristics [74]. However, the mechanism of these three genes located near SSC12:50896936 on porcine BFT traits needs to be further studied.
Candidate genes for CLP and CFP traits
For the traits of CLP and CFP, three loci, including SSC10:46524739, SSC11:47266057 and SSC11:47266119 were significantly associated with CLP and CFP traits. SSC10:46524739 could explain 8.98%, and 10.873% of the PV for CLP and CFP, respectively. The SNP locus was located in the intron region of FAM171A1. GO analysis showed that the FAM171A1 gene was significantly enriched in regulation of cell shape (Table S3). The FAM171A1 gene is a member of the family of sequence similarities, which is responsible for encoding the APCN/FAM171A1 protein. APCN/FAM171A1 is an evolutionarily conserved 98 kDa transmembrane type I glycoprotein expressed in various cells, which is involved in the regulation of the cytoskeletal dynamics and thereby the cell shape [75]. According to its function and extensive expression in various cells, perhaps FAM171A1 could be considered a candidate gene for the porcine CLP and CFP traits and further studied its effects on the growth, proliferation and differentiation of muscle cells and fat cells. However, no genes were found within 100kb upstream and downstream of the SSC11:47266057 and SSC11:47266119, but these two loci could explain about 10% of the PV for CLP and CFP, respectively. It was inferred that these two SNPs might be important for the CLP and CFP of the pigs and might be potential targets for further molecular breeding research.
Six genes, including SDS, LHX5, SDSL, PLBD2, DTX1, and RASAL1 were found nearby SSC14:38533013 associated with CLP. Among these six genes, the intron region of PLBD2 contained the SNP loci. GO analysis showed that PLBD2 was enriched in lipid catabolic process and hydrolase activity. Wang et al. find that PLBD2 and other nine genes, predominantly participate in the lipid catabolic process, which is in agreement with the favoured selection for fatness in Chinese pigs through analysis of the selected regions for Tongcheng pigs using the selective sweep method [76]. The PLBD2 gene should be considered a strong candidate gene for porcine CLP traits. Finally, the SNP of SSC10:50888840 related to CFP traits was located in the intron region of KIAA1217. A study finds that the KIAA1217 gene shows differences in promoter methylation and mRNA expression in the omental visceral adipose tissue between non-obese and obese individuals [77]. Perhaps, the KIAA1217 gene might be used as a candidate gene for the porcine CFP trait to further study its effects on lipid metabolism.
However, these identified loci and genes need to be further verified in more pig populations, and their functions also need to be validated by more biological experiments in pigs.