==== Front PLoS Genet PLoS Genet plos PLOS Genetics 1553-7390 1553-7404 Public Library of Science San Francisco, CA USA 37339141 10.1371/journal.pgen.1010820 PGENETICS-D-22-01410 Research Article Biology and Life Sciences Genetics Single Nucleotide Polymorphisms Biology and Life Sciences Organisms Eukaryota Animals Vertebrates Amniotes Mammals Swine Biology and Life Sciences Zoology Animals Vertebrates Amniotes Mammals Swine Biology and Life Sciences Computational Biology Genome Analysis Genome-Wide Association Studies Biology and Life Sciences Genetics Genomics Genome Analysis Genome-Wide Association Studies Biology and Life Sciences Genetics Human Genetics Genome-Wide Association Studies Biology and Life Sciences Genetics Genomics Biology and Life Sciences Genetics Biology and Life Sciences Genetics Genetic Loci Quantitative Trait Loci Biology and Life Sciences Molecular Biology Molecular Biology Techniques Gene Mapping Restriction Fragment Mapping Electrophoretic Mobility Shift Assay Research and Analysis Methods Molecular Biology Techniques Gene Mapping Restriction Fragment Mapping Electrophoretic Mobility Shift Assay Biology and Life Sciences Computational Biology Genome Analysis Biology and Life Sciences Genetics Genomics Genome Analysis Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs Integrated analysis of gwas and 3D genomics reveal BMP2 regulating LMD in pigs Miao Yuanxin Conceptualization Formal analysis Funding acquisition Software Writing – original draft 1 2 Zhao Yunxia Conceptualization Formal analysis Methodology Writing – original draft 1 Wan Siqi Formal analysis Software Writing – original draft 1 Mei Quanshun Software 1 Wang Heng Writing – review & editing 1 Fu Chuanke Software 1 Li Xinyun Writing – review & editing 1 3 Zhao Shuhong Funding acquisition Writing – review & editing 1 3 Xu Xuewen Methodology Writing – review & editing 1 * https://orcid.org/0000-0001-6134-2627 Xiang Tao Conceptualization Funding acquisition Methodology Writing – original draft 1 * 1 Key Laboratory of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education & Key Laboratory of Swine Genetics and Breeding of Ministry of Agriculture, Huazhong Agricultural University, Wuhan 430070, China 2 Research Institute of Agricultural Biotechnology, Jingchu University of Technology, Jingmen 448000, China 3 Hubei Hongshan laboratory, Huazhong Agricultural University, Wuhan 430070, China Groenen Martien Editor Wageningen University & Research, NETHERLANDS The authors have declared that no competing interests exist. * E-mail: Xuewen_Xu@mail.hzau.edu.cn (XX); Tao.Xiang@mail.hzau.edu.cn (TX) 20 6 2023 6 2023 19 6 e101082010 12 2022 7 6 2023 © 2023 Miao et al 2023 Miao et al https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Background The lack of integrated analysis of genome-wide association studies (GWAS) and 3D epigenomics restricts a deep understanding of the genetic mechanisms of meat-related traits. With the application of techniques as ChIP-seq and Hi-C, the annotations of cis-regulatory elements in the pig genome have been established, which offers a new opportunity to elucidate the genetic mechanisms and identify major genetic variants and candidate genes that are significantly associated with important economic traits. Among these traits, loin muscle depth (LMD) is an important one as it impacts the lean meat content. In this study, we integrated cis-regulatory elements and genome-wide association studies (GWAS) to identify candidate genes and genetic variants regulating LMD. Results Five single nucleotide polymorphisms (SNPs) located on porcine chromosome 17 were significantly associated with LMD in Yorkshire pigs. A 10 kb quantitative trait locus (QTL) was identified as a candidate functional genomic region through the integration of linkage disequilibrium and linkage analysis (LDLA) and high-throughput chromosome conformation capture (Hi-C) analysis. The BMP2 gene was identified as a candidate gene for LMD based on the integrated results of GWAS, Hi-C meta-analysis, and cis-regulatory element data. The identified QTL region was further verified through target region sequencing. Furthermore, through using dual-luciferase assays and electrophoretic mobility shift assays (EMSA), two SNPs, including SNP rs321846600, located in the enhancer region, and SNP rs1111440035, located in the promoter region, were identified as candidate SNPs that may be functionally related to the LMD. Conclusions Based on the results of GWAS, Hi-C, and cis-regulatory elements, the BMP2 gene was identified as an important candidate gene regulating variation in LMD. The SNPs rs321846600 and rs1111440035 were identified as candidate SNPs that are functionally related to the LMD of Yorkshire pigs. Our results shed light on the advantages of integrating GWAS with 3D epigenomics in identifying candidate genes for quantitative traits. This study is a pioneering work for the identification of candidate genes and related genetic variants regulating one key production trait (LMD) in pigs by integrating genome-wide association studies and 3D epigenomics. Author summary The loin muscle depth (LMD) is positively correlated with lean meat content and can therefore be used to estimate the lean meat content of pigs. The loin muscle depth is a quantitative trait that is influenced by whole-genome distributed polygenes. To identify the possible candidate genes and mutations impacting the LMD, we used genome-wide association study (GWAS), which is commonly used to identify candidate genomic regions associated with complex growth traits, to mapping the genomic regions. As over 90% of the identified variants have been localized to non-coding regions of the pig genome and their functions in phenotype regulation are poorly understood, we further combined 3D epigenomics to understand the mechanism of genetic variants regulating LMD. In combination with different types of data and analyses, we identified a 10 kb quantitative trait locus (QTL) on chromosome 17 that was significantly associated with LMD. We identified the BMP2 gene as a major candidate gene, and SNP rs1111440035 and rs321846600 were identified as likely candidate mutations affecting the LMD. Our study is unique in its attempt to identify candidate genetic variants by integrating GWAS and 3D epigenomics in pigs. the National Key Research and Development Program of China 2021YFF1000601 & 2022YFD1301900 https://orcid.org/0000-0001-6134-2627 Xiang Tao http://dx.doi.org/10.13039/501100018531 Major Science and Technology Projects in Yunnan Province No.2020ABA016 https://orcid.org/0000-0001-6134-2627 Xiang Tao the China Agriculture Research System of MOF and MARA CARS-35 Zhao Shuhong Natural Science Foundation in Hubei Province No. 2022CFB021 Miao Yuanxin the Scientific Research Program Key Project of Hubei Provincial Department of education D20214301 Miao Yuanxin the Natural Science Foundation of Jingmen City 2022YFZD051 Miao Yuanxin the research fund from Jingchu University of Technology ZD202102, 2022ZD005 Miao Yuanxin TX acknowledges funding from the National Key Research and Development Program of China (No. 2021YFF1000601 & 2022YFD1301900), Major Science and Technology Projects in Hubei Province (No.2020ABA016); SZ acknowledges funding from the China Agriculture Research System of MOF and MARA (CARS-35); YM acknowledges funding from Natural Science Foundation in Hubei Province (No. 2022CFB021), the Scientific Research Program Key Project of Hubei Provincial Department of education (D20214301), the Natural Science Foundation of Jingmen City (2022YFZD051), and the research fund from Jingchu University of Technology (ZD202102, 2022ZD005). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. PLOS Publication Stagevor-update-to-uncorrected-proof Publication Update2023-06-30 Data AvailabilityThe data of phenotype and genotype are deposited in supporting information files S1 Data and S2 Data. Data Availability The data of phenotype and genotype are deposited in supporting information files S1 Data and S2 Data. ==== Body pmcIntroduction Lean pigs provide the majority of pork for the consumer market, and lean meat content is a critical breeding goal for the pig industry [1,2]. The loin muscle depth (LMD) is defined as the minimum distance from the vertebral channel to the cranial end of the gluteus medius muscle. The LMD is usually measured by ultrasound, and it is significantly positively correlated with lean meat content. Therefore, selection for LMD enhances the lean meat content of pigs over time. Revealing the genetic mechanisms governing LMD and identifying genes significantly associated with LMD may further enhance the efficiency of LMD improvement. Genome-wide association study (GWAS) is a powerful and effective strategy for detecting QTL regions associated with complex traits [2–6]. It has been successfully used to identify candidate genes associated with both LMD and loin muscle area (LMA). LMD is significantly positively correlated with LMA (r = 0.945), and it also can be used to predict lean meat content [7]. For instance, 9 independent genome-wide SNPs accounting for 27.51% of phenotypic variation associated with LMD were identified from 370 Chuying-black pigs using GWAS [8]. Another example has demonstrated that GWAS meta-analysis has identified multiple QTL regions related to LMA or LMD across 6,043 Duroc pigs, and these QTL regions contained 75 significantly associated SNPs [9]. To date, 415 QTLs for loin muscle area and 97 QTLs for loin muscle depth are available in the Pig QTL Database (https://www.animalgenome.org/cgi-bin/QTLdb/SS/index) [6,10–12]. Although over 90% of these identified variants have been mapped to non-coding regions in the pig genome, there is limited understanding of their functions in these non-coding regions [13]. To uncover the molecular mechanisms of these variants, it is necessary to elucidate not only the genetic variations but also the epigenetic activities accompanying this genetic variation. Previous studies have indicated that integrating GWAS and epigenetic data, including chromatin states [14], open chromatin regions [15], high-throughput chromosome conformation capture (Hi-C) [16], and RNA-seq [17], could enhance our understanding of the identified variants and their target genes as well as the relationships between genotypes and phenotypes. Recently, the accumulation of epigenetic data [18] offers an opportunity to interpret GWAS results from pigs to reveal pertinent genetic mechanisms of specific traits. This study was aimed at identifying major candidate genes as well as key candidate genetic variants influencing phenotypic variation in LMD and exploring the genetic mechanisms governing LMD regulation in Yorkshire pigs. In this study, we performed a GWAS of LMD in Yorkshire pigs, followed by a chromatin state analysis in a region significantly associated with LMD based on genome-wide high-resolution profiles from ChIP-seq and Hi-C data. Our work provides a foundational framework and a successful case for combining GWAS and 3D epigenomics in pigs. Materials and methods Ethics statement All procedures involving tissue samples collection and animal care were performed according to the approved protocols and ARRIVE (Animal Research: Reporting In Vivo Experiments) guidelines and were approved by the Ethics Committee of Huazhong Agricultural University (HZAUSW-2018-008). Phenotypic recordings and genotypes In this study, loin muscle depth (LMD) from the 10th rib to 11th rib of pigs was measured in Yorkshire pigs with the weight of 100 ± 5 kg by an Aloka 500V SSD B ultrasound. In total, LMD of 16,533 pigs was recorded in the years from March 2014 to January 2018. The pedigrees of these measured pigs could be traced back at least five generations with a total of 42,245 pigs included in these pedigrees. Genomic selection was started late in 2016 and finished at the end of 2018 due to the outbreak of African swine fever. Ear tissues were collected by the criterion of at least 2 males and 2 females in each litter. A total of 2,733 Yorkshire (YY) pigs provided both LMD recordings and ear tissue samples simultaneously. Total DNA was extracted from these 2,733 pigs and genotyped using the Geneseek Porcine 50K SNP Chip (Neogen, Lincoln, NE, United States), and 50,703 SNP markers across the genome were obtained. SNPs were mapped to the Sscrofa11.1 pig genome assembly. Quality control (QC) was carried out with the following parameters set: individual call rate≥90%; SNP call rate≥90%; minor allele frequency≥0.01; Hardy-Weinberg equilibrium P-value≥10−6. After quality control, a total of 42,772 SNPs were obtained from each of 2,733 pigs, and these SNPs were used for subsequent genome-wide association analysis. Pre-correction of phenotypes The loin muscle depth (LMD) of 16,533 pigs was pre-corrected. The single-step GBLUP (ssGBLUP) method, which incorporates both pedigree and genomic information, was used to estimate variance components and genomic estimated breeding values (GEBVs) [19,20]: y=Xb+Za+e (1) where y is a vector of original phenotypic values of LMD of 16,533 pigs; Xb accounts for the fixed effects of herd-year-season, sex, and age; a is a vector of random additive genetic effects; and Z is the incidence matrix; e is a vector of residual effects, which was assumed to follow normal distribution as e∼N(0,Iσe2), where I is an identity matrix. It was assumed that the random additive effects follow a normal distribution a∼N(0,Hσa2), where H is a relationship matrix combining pedigree and genomic information, which was constructed by previously reported method [19]. σa2 is the additive genetic variance. Software DMU [21] was used to estimate the variance components and solve the genomic model. The corrected phenotypes yc were calculated as the sum of estimated additive genetic effects and the residuals (yc=a^+e^) for all the pigs. Genome-wide association studies The genome-wide association study (GWAS) was performed on 2,733 pigs with both phenotypic and genotypic information by using a mixed linear model-based association analysis (MLMA) in software rMVP [22]. The following mixed linear model was used for GWAS: yc=1μ+Xg+Wu+e, (2) where yc is a vector of corrected phenotype of LMD in the genotyped 2,733 Yorkshire pigs; μ is the overall mean of corrected phenotypes; 1 is a vector of ones; X is a matrix of the SNP genotypes with entry 0, 1, 2 indicating genotype AA, AB, and BB, respectively; g is the fixed additive genetic effect of the analyzed SNP; u is a vector of random polygenic effects, with an assumption that u∼N(0,Gσu2), where G is the marker-based additive genomic relationship matrix, which was constructed by the method reported by Vanraden (2008) [23]; σu2 is the polygenic variance; W is the incidence matrix between the corrected phenotype and the corresponding random polygenic effects; e is a vector of random residual effects, with an assumption that e∼N(0,Iσe2). In Bonferroni corrections, the genome-wide significant threshold was set as (−log10[0.05/number of SNPs] = 5.93). Identification of LD block and QTL analysis Identification of linkage disequilibrium (LD) blocks was performed in the chromosomal regions containing the identified significantly associated SNPs by software Haploview [24]. The QTLs located in the identified LD blocks were searched from Pig QTL database (pigQTLdb, https://www.animalgenome.org/cgi-bin/QTLdb/SS/index). Haplotype analysis Haplotypes were constructed in the LMD-associated region on SSC17 through the Linkage disequilibrium and linkage analysis (LDLA) method [25]. The LDLA method can accurately localize the identified QTL regions by combining the results of populational linkage disequilibrium and within-family linkage, thus LDLA was used for the haplotype analysis [26,27]. The haplotypes were constructed by the software PHASEBOOK [28] using the Hidden Markov Model [28] based on an assumption that there existed a predetermined number of ancestral haplotypes (K = 20), and that all haplotypes in the population were derived from these ancestral haplotypes [29]. A likelihood ratio test (LRT) was performed along chromosomes to test the presence of a QTL region in the given map [30]. The 95% confidence interval (CI) was calculated as the LRT value of the most significant loci minus 2 [31,32]. Integrated analysis of pig epigenomics datasets The topologically associated domains (TAD), sub-domains, Hi-C contact matrix data, significant H3K27ac peaks, and enhancer-gene pairs were all from our previous study based on the Sscrofa11.1 pig genome [18]. Hi-C data were analyzed using HiC-Pro pipeline version 2.9.0 software for genome mapping. The insulation score and top domain methods were used to perform TAD calling and sub-domain identification. ChIP-seq data analysis included reads mapping (BWA v0.7.15), low-quality reads filtering (SAMTools v1.9 and Picard v1.126), and peak calling (MACS2 v2.1.0). The details of Hi-C and ChIP-seq data analysis were described in our previous study [18]. In this study, the Hi-C contact and TAD structure integrated with GWAS were visualized using Juicebox [33]. The significant H3K27ac peaks across various tissues were merged using BEDTools v2.26.0 [34]. Significant H3K27ac peaks surrounding the candidate genes were visualized in the IGV browser [35]. SNP polymorphism identification and association analysis To further identify important functional variants, the genomic region between 15.51 Mb and 16.31 Mb on SSC17 was deeply sequenced (>20X) by using target region sequencing technology on 732 randomly selected individuals across the entire population. Target region sequencing libraries were created from isolated DNA according to the manufacturer’s instructions. The high-quality libraries were sequenced using the Illumina HiSeq3000 platform, which generated paired-end sequencing data. Quality control of paired-end sequencing reads was conducted with Trimmomatic, followed by alignment to the pig reference genome using BWA. The variants were identified using GATK software according to specific criteria (Qual score ≥ 30, QD < 20.0, ReadPosRankSum < −8.0, FS > 10.0 and QUAL < $MEANQUA). Subsequently, SNPs located between 15.51 Mb and 16.31 Mb on SSC17 and distributed within the region of significant H3K27ac peaks were selected for trait association analysis with the corrected LMD phenotypes across the 732 randomly-selected individuals. The association analysis between an individual SNP and LMD was carried out using a generalized linear model using R software. This model was Y = 1μ+genotype+e, where Y is the response vector of the LMD, μ is the mean of LMD, genotype contains three different levels of genotypes (AA, AB, and BB). The genotype was considered a fixed effect in the model, and e represents a vector of residual errors. For LMD, the least square means of different genotypes (AA, AB, and BB) were compared through Least Significant Difference (LSD) in R. An SNP with a p -value< 0.05 was considered to be significantly associated with LMD. Dual-luciferase expression assays of promoter and enhancer regions For the SNPs that were identified to be significantly associated with LMD through the aforementioned steps, JASPAR [36] was used to identify transcription factor binding sites. For the five SNPs found in the enhancer region of BMP2, which impact transcription factor binding sites (S1 Fig), we conducted dual-luciferase expression assays to verify their functions. In the enhancer region, the 400-bp genome region flanking the significantly associated SNPs was cloned. As the distance between the SNP rs328487632 and rs319025934 is 24bp, we only cloned a single 800bp genome region containing these two SNPs. In total, four fragments in the enhancer region were cloned. Dual-luciferase expression experiments were performed for the promoter region of candidate genes, where the region (2.5 kb upstream) surrounding the transcription start site (TSS) was considered the promoter of the candidate genes. Genomic DNA was extracted from the ears of Yorkshire pigs. The promoter and enhancer regions of candidate genes were amplified from genomic DNA through PCR using Phanta Max Super-Fidelity DNA Polymerase (P505, Vazyme). The PCR product was subsequently cloned into the pGL3- basic (E1751, Promega) and pGL3-promoter vectors (E1761, Promega) upstream of the luciferase gene for promoter and enhancer assays, respectively, using restriction enzymes NcoI (FD0574, Thermo) and KpnI (FD0524, Thermo). Dual-luciferase assays were conducted in PK15 (Porcine Kidney 15) cells which were cultured in a 37°C incubator with 5% CO2. PK15 cells were plated into 48-well plates and co-transfected with reporter vectors and pRL-TK renilla luciferase control vector (E2241, Promega) using Lipofectamine 2000 (11668500, Invitrogen). The transfected cells were incubated for 24 hours prior to lysis, which were lysed by using the Dual-Luciferase Reporter Assay System (11402ES60, Yeasen). For each construct, the reporter assay was repeated at least three times independently, and the results from a single representative experiment are presented in this study. Statistical analysis was conducted using a Student’s t-test, and a p-value lower than 0.05 was considered significant. Electrophoretic mobility shift assays (EMSA) PK15 nuclear extracts were prepared using NE-PER Nuclear and Cytoplasmic Extraction Reagents (78833, Thermo Fisher Scientific) and quantified using the BCA method (P0010S, Beyotime). Single-stranded and reverse-complement DNA probes were synthesized with or without a 5′-end biotin label. For the candidate SNPs, probes were designed as follows: for SNP rs1111440035 (M4 in the promoter, chr 17:15749990 bp), the M4-L-CC probe forward strand was 5′-CCCACCCGAACGACCTCGGGGCGA-3′; M4-H-GG probe forward strand was 5′-CCCACCCGAAGGACCTCGGGGCGA-3′; for SNP rs80791204 (M5 in the promoter, chr 17:15750750 bp), the M5-H-GG probe forward strand was 5′-AGGGAGAATAACTTGGGCTCCTCACTTCGCG-3′; and the M5-L-CC probe forward strand was 5′-AGGGAGAATAACTTGCGCTCCTCACTTCGCG-3′; for SNP rs321846600 (in the enhancer), the TT probe forward strand was 5′-GACAACCAGATCCATCTGGGCACCAGTC-3′; and the CC probe forward strand was 5′-GACAACCAGATCCACCTGGGCACCAGTC-3′. An EMSA assay was performed using a Light Shift Chemiluminescent EMSA Kit (89880, 20148E, Thermo Fisher Scientific) according to the product manuals. We prepared and performed the binding reactions as follows: for the negative groups (lanes 1 and 4), 1 μL of 10× binding buffer, 0.5 μL of 1 μg/ μL poly(dI-dC), and 6.5 μL of ddH2O; for the experiment groups (lanes 2 and 5), 1 μL (2 μg) of nuclear protein, 1 μL of 10× binding buffer, 0.5 μL of 1 μg/ μL poly(dI-dC), and 6.5 μL of ddH2O; for the competition groups (lanes 3, 6, 7, and 8), 1 μL (2 μg) of nuclear protein, 1 μL of 10× binding buffer, 0.5 μL of 1 μg/ μL poly(dI-dC), 2 μL of 1 pmol/μL unbiotin-labeled probes, and 4.5 μL of ddH2O. Binding reactions were incubated at room temperature for 20 minutes, followed by an additional incubation at room temperature for 30 minutes after the addition of biotin-labeled probes. Then, 2.5 μL of 5× Loading Buffer was added to each 10μL reaction. To run the 6% polyacrylamide gel, the voltage was set to 100 V, and samples were electrophoresed until the bromophenol blue dye migrated approximately 3/4 through the length of the gel in 0.5× TBE buffer. We transferred the bands at 380 mA (~100V) for 60 minutes. Finally, the biotin-labeled DNA was detected using chemiluminescence. Results Descriptive statistical analysis and genetic parameters In this study, we obtained 16,533 recordings of LMD in Yorkshire pigs. Descriptive statistics of the LMD recordings are listed in Table 1. The mean of LMD was 6.19 cm with a standard deviation (SD) of 0.74 cm. The additive genetic variance of LMD was 0.095, with a standard error of 0.016. The residual variance was 0.160, with a standard error of 0.035. The estimated heritability of LMD was 0.373, with a standard error of 0.018. 10.1371/journal.pgen.1010820.t001 Table 1 Descriptive statistics of loin muscle depth in the Yorkshire population. Traits #indiv Min Mean Max SD Loin muscle depth (cm) 16533 1.52 6.19 8.94 0.74 Five significant SNPs associated with LMD were identified by GWAS A linear mixed model analysis was applied to perform the GWAS analysis of the LMD trait. The significance threshold was calculated as the cut-off after the Bonferroni correction. In total, 5 SNPs reached the significance threshold of 5.93 (−log10(0.05/42772) = 5.93) (Fig 1). All significantly associated SNPs were located within the region of 15.51 to 16.31 Mb on SSC17 (with a span of 0.8 Mb). Detailed information regarding the SNPs and their nearest genes is shown in Table 2. 10.1371/journal.pgen.1010820.g001 Fig 1 (a) Manhattan plot and (b) QQ plot showing genome-wide association for loin muscle depth (LMD) in Yorkshire pig. The red dotted line in (a) represents -log10(p)-value = 5.9. 10.1371/journal.pgen.1010820.t002 Table 2 Summary of significantly associated SNPs for loin muscle depth (LMD) and relevant genes in Yorkshire pigs. SNP SSC Position(bp) Alleles P value Region Nearest Gene WU_10.2_17_16896163 17 15518001 C/T 1.29E-08 Intergenic ENSSSCG00000043546 WU_10.2_17_16951872 17 15534531 G/A 1.21E-08 Intergenic ENSSSCG00000043546 WU_10.2_17_17013787 17 15710331 C/T 5.31E-07 Intergenic BMP2 MARC0112426 17 15755711 C/A 4.47E-12 intronic BMP2 WU_10.2_17_18106530 17 16312580 T/C 2.76E-08 Intergenic ENSSSCG00000025527 LMD-related SNPs are located in the same TAD Our previous study identified the boundaries of topologically associated domains (TADs) using Hi-C data from pig muscle tissues. These boundaries served to limit chromatin interactions, mediated by specific proteins, between different TADs in nuclei. The interaction effects of cis-regulatory elements within the same TADs were significantly stronger than those spanning the two nearest adjacent TADs [18]. Based on this finding, we integrated the GWAS results and Hi-C data to identify candidate interaction regions surrounding the SNPs that were significantly associated with LMD. Results suggested that these significantly associated SNPs (spanning 0.8 Mb and ranging from 15.51 Mb to 16.31 Mb on SSC17) were all located in the same TAD region (SSC17: 15.08 to 16.76 Mb). Three sub-domains covered all the significantly LMD-associated SNPs. Among these three sub-domains, one sub-domain (15.65 to 15.89 Mb) was completely within the 0.8-Mb significantly LMD-associated SNP region (15.51 to 16.31 Mb), while the other two were partially within it (Fig 2B). These results demonstrated that the significant SNP-located region had limited opportunity to interact with genome regions located in other TADs due to the presence of boundaries between different TADs. 10.1371/journal.pgen.1010820.g002 Fig 2 Association results for SSC17. (a) Manhattan plot of genome-wide p-values for pig chromosome 17. (b) Hi-C contact heatmap surrounding a significantly associated QTL region at 10 Kb resolution. (c) Map of characterized genes from 15.51 to 16.31 Mb on SSC17. We further investigated genes located in the aforementioned TAD region (15.08–16.76 Mb) on pig chromosome 17 and uncovered 11 genes located in this region. Of these 11 genes, seven were long non-coding RNA genes, whereas the other four genes were coding genes. Detailed information is shown in Fig 2 and S1 Table. This TAD region was investigated in the pigQTLdb database, showing that QTLs in this genomic region (15.08–16.76 Mb on SSC17) were associated with body weight, body height, intramuscular fat content, average daily gain, average backfat thickness, carcass length, and other traits (S2 Table), which may influence the loin muscle depth. Linkage disequilibrium analysis highlights the candidate QTL region harboring BMP2 To further characterize the candidate QTL region within the TAD (SSC17: 15080000–16760000), we investigated linkage disequilibrium (LD) patterns around the significantly LMD-associated SNPs. Four LD blocks were detected in this region using the confidence interval algorithm in Haploview software (Fig 3A). MARC0112426, the most significantly associated SNP, was located in the intron of BMP2 (SSC17: 15750487–15762982) within LD block 2 (15.65 Mb to 15.75 Mb). Two significantly associated SNPs (WU_10.2_17_17013787 and MARC0112426) were completely linked (D’ = 1) in LD block 2. We further studied the Hi-C contact map surrounding LD block 2. There was a complete sub-domain (15.65 to 15.89 Mb) covering LD block 2 (15.65 to 15.75 Mb) (Fig 3B), which indicated LD block 2 has a high opportunity to interact with genomic regions within this sub-domain. Although the SNPs located inside and outside of LD block 2 showed a high linkage disequilibrium (LD) with each other, the sub-domains isolate this disequilibrium. Thus, SNPs outside the sub-domains were considered to be irrelevant to the LMD trait. Overall, these results indicated that LD block 2 is a compelling candidate QTL region for the LMD trait. 10.1371/journal.pgen.1010820.g003 Fig 3 Linkage analysis (LDLA) integrating Hi-C data locating target regions significantly associated with LMD. (a) Linkage disequilibrium blocks were detected in the region (15.51 to 16.31 Mb on SSC17). SNPs in red boxes have the highest P value. As the chip is named in accordance with version 10.2, the figure uses version 10.2. The location has been converted to version 11.1 in subsequent analyses. (b) Hi-C contact heatmap of significant SNP locations (LMD significantly associated QTL region 15.51–16.31 Mb). (c) Likelihood ratio test (LRT) profiling using the combined linkage disequilibrium and linkage analysis (LDLA) approach for LMD. The x-axis represents the scan steps of the analysis, and the y-axis represents the LRT value. The top six LRT values are denoted in red, and the 95% confidence interval (CI) defined by the LRT-dropoff-2 method is indicated in a red-shaded block. Furthermore, LDLA results demonstrated that the block 2 region (15.65 to 15.75 Mb) was significantly associated with LMD (Fig 3C and S3 Table). The SNPs with the top six LRT values were within the 95% confidence interval of the highest LRT values [31,32], and the region (15.65 to 15.89 Mb) containing these six SNPs completely covered the identified LD block2 region (15.65 to 15.75 Mb). The SNP with the highest LRT value (WU_10.2_17_17075196) was located at SSC17:15689085. These results confirmed that the LD block2 region (SSC17:15659761–15755711) was an important candidate QTL region associated with LMD, aligned with the GWAS and Hi-C integrated results. Only one gene, BMP2, was located in this candidate QTL region, suggesting that BMP2 is a candidate gene related to LMD (Fig 3C). SNPs in cis-regulatory elements of BMP2 are significantly associated with LMD Cis-regulatory elements of the pig genome were identified in our previous study [18]. To further comprehend the genetic mechanisms of LMD, we investigated the association between SNPs in the cis-regulatory elements of the pig genome and the LMD trait (Fig 4). 10.1371/journal.pgen.1010820.g004 Fig 4 Visualization of SNPs using different methods or validation locating in significant H3K27ac ChIP-seq peak regions from various tissues [18]. (a) Hi-C contact heatmap in the candidate QTL region (LD block 2, chr 17: 15659761–1575571) of LMD. (b) Significant enhancer-gene pairs in the candidate QTL region of LMD [18]. (c) SNPs from the Yorkshire pig whole genome sequence data, and target sequencing validated in H3K27ac ChIP-seq peak regions of various tissues. The purple box and sky blue box represent the significant (P<0.05) SNPs validated by target sequencing. (d) Scatter plot showing single-SNP trait association analysis of targeted sequencing data. The red dotted line in (d) represents a p-value = 0.05. To further identify the functional variants, we analyzed SNPs using target region sequencing technology. For the 732 randomly selected individuals, 4,980 SNPs were captured in the GWAS candidate region (15.51 Mb to 16.31 Mb on SSC17). According to the association analysis, there were 184 SNPs significantly associated with LMD (Fig 4D and S4 Table), of which five SNPs, chr17:15674366 (rs321846600), chr17:15683553 (rs328487632), chr17:15683577 (rs319025934), chr17:15684170, and chr17:15724776 (rs321766789), were located in cis-regulatory elements of the pig genome. These five SNPs were distributed in the enhancer region of BMP2. Moreover, the Hi-C contact map and enhancer-gene correlation results suggested that enhancers in the LMD QTL regions interacted with the BMP2 promoter (Fig 4A and 4B). These results suggested that the identified SNPs in the cis-regulatory elements of BMP2 are associated with LMD in pigs. Important candidate variant scanning in the enhancer region of BMP2 Based on the target region sequencing results, five significantly associated SNPs were located in the BMP2 enhancer region. SNPs in regulatory regions may function by modulating transcription factor binding. Therefore, motif analyses were performed to further validate the functional consequences of these five SNPs. Results suggested that these five SNPs may all disturb the TF binding motif (S1 Fig). We hypothesized that the genomic regions harboring these five SNPs are able to regulate the LMD by altering enhancer function. To validate our hypothesis, the 400-bp genomic region flanking each significantly associated SNPs with different allelic combinations was cloned into pGL3-promoter luciferase reporter vectors respectively. Because the distance between SNP rs328487632 and rs319025934 is only 24 bp, we cloned one 800-bp genome region containing these two SNPs. From reporter assays, we observed stronger luciferase activity for each cloned genomic region compared to the empty vector, supporting the hypothesis that these five SNPs can function as enhancers (Fig 5A–5D). Furthermore, for the cloned 800-bp genome region surrounding the five aforementioned SNPs, we compared the luciferase activities between alleles. Results showed that the enhancer activity for SNPs rs328487632&rs319025934, and rs321766789 did not show a significant change between alleles (Fig 5A and 5B). The cloned 800-bp genomic region covering SNP 15684170 contains more than one SNP in the region. To exclude effects on enhancer activity from other SNPs, in addition to comparing the enhancer activity between the T and G substitution at position 15684170, we mutated the T to a G in position 15684170. Results (Fig 5D) showed that the enhancer activity did not change significantly, indicating that SNP 15684170 is not an important site for enhancer activity. 10.1371/journal.pgen.1010820.g005 Fig 5 Dual-luciferase reporter and EMSA assays. (a–d) Luciferase reporter assays containing the SNPs 15683553(rs328487632), 15683577(rs319025934), 15724776(rs321766789), 15684170, and 15674366(rs32184660), for the enhancer assay. (d) 15684170T-G represents the SNP 15684170 mutating a T into a G. Expression levels were measured after transfection into the porcine kidney cell line, PK15. The pGL3 promoter vector containing the SV40 promoter was used for the enhancer assay. Luciferase signals were normalized to Renilla signals (n = 3). Data are presented as mean ± SEM, and p values were calculated using a Student’s t-test. * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001. (e) EMSA results of rs321846600. U-TT represents the unlabeled TT probe, similar to U-CC, and B-TT represents the biotin-labeled TT probe, similar to B-CC. Three replicates were performed for each experiment. Among these five SNPs, the presence of a C at rs321846600 showed significantly higher enhancer activity than a T (Fig 5C). This implied that the SNP rs321846600 in the enhancer region of BMP2 may be important for the regulation of BMP2 expression. Thus, we considered the SNP rs321846600 as a candidate variant that may be functionally related to the LMD trait. Variant scanning in the promoter of BMP2 identified two functional mutations We further characterized the BMP2 promoter and investigated its polymorphisms. We cloned the promoter region of porcine BMP2 from 10 Yorkshire pigs to scan for potential variants, and five variants (4 SNPs and one Indel) were identified. These five variants are completely linked and formed two haplotypes: T(CAAAC)TGG (BMP2 H) and C(T----)ACC (BMP2 L) (Fig 6A). The promoter with the haplotype T(CAAAC)TGG (BMP2 H) had higher transcription activity than C(T—-)ACC (BMP2 L), as found by the dual-luciferase reporter assays (Fig 6B). 10.1371/journal.pgen.1010820.g006 Fig 6 Identification of functional mutations in the porcine BMP2 promoter. Luciferase reporter assays using vectors containing the SNP locus for the promoter assay. (a)(c)(e) SNP location of each vector. The transcriptional start site (TSS) was defined as +1. M1 represents the vector containing mutations in the first sites. (b)(d)(f) Luciferase reporter assays using vectors containing the SNP locus for promoter assay. The pGL3 basic vector lacking a promoter was used for the promoter assay. Luciferase signals were normalized to Renilla signals. Each assay involved three independent experiments. The results from one representative experiment in the porcine kidney cell line, PK15, are shown. Data are presented as mean ± SEM, and p values were calculated using a Student’s t-test. * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001. To identify important candidate variants influencing promoter activities between BMP2 H and BMP2 L, we created mutant reporter vectors (BMP2 H>L and BMP2 L>H) based on the five identified SNPs (Fig 6C). The reporter assay suggested that the promoter activity had changed significantly (Fig 6D), which suggested that M4 (rs1111440035, chr 17:15749990 bp) or M5 (rs80791204, chr 17:15750750 bp) altered the promoter activities of BMP2. The original alleles in the M4 (BMP2 M4 H>L) and M5 (BMP2 M5 H>L) of BMP2 H were mutated to the corresponding alternative alleles one nucleotide at a time (Fig 6E). Results showed that when the alleles in either the M4 (G>C) or M5 (G>C) sites of BMP2 H were mutated, the promoter activity significantly decreased (Fig 6F). Therefore, we concluded that both rs1111440035 and rs80791204 caused different promoter activities of BMP2 H and BMP2 L. To further validate the identified candidate SNPs, motif analysis was conducted for the two SNPs. Motif analysis showed that rs1111440035 (M4, C>G) disturbed the binding of the transcription factor NR2C1 (Fig 7A). Similarly, rs80791204 (M5, G>C) overlaps with a predicted HLTF binding motif (Fig 7B). Moreover, multiple sequence alignments (Fig 7C and 7D) revealed that the rs1111440035(M4) within the porcine BMP2 promoter was conserved across mammals. Genomic Evolutionary Rate Profiling (GERP) scores are often used to measure the conserved nature of gene sequences across species, and a high GERP score suggests that sequence is highly conserved. We added GERP scores to determine the conservation of M4 and M5. For M4, the GERP score is 1.91; for M5, the GERP score is −5.13. These GERP scores shows that M5 is less conserved than M4. 10.1371/journal.pgen.1010820.g007 Fig 7 Analysis and EMSA assay for rs1111440035 and rs80791204. (a)(b) Motif analysis for rs1111440035(M4) and rs80791204(M5). (c)(d) Multiple sequence alignment across mammalian DNA elements surrounding the M4 and M5 mutation sites. (e)(f) The EMSA results. Here, “protein” represents PK15 nuclear extracts. U-M4-L-CC represents unlabeled M4-L-CC probe, similar to U-M4-H-GG, U-M5-L-CC, and U-M5-H-GG, while B-M4-L-CC represents biotin-labeled M4-L-CC probe, similar to B-M4-H-GG, B-M5-L-CC, and B-M5-H-GG. Three replicates were performed for each experiment. Verification of important candidate variants by EMSA Next, we performed EMSAs using nuclear proteins from PK15 cells to evaluate binding affinities for important candidate variants (SNP rs321846600 in the enhancer region; rs1111440035 and rs80791204 in the promoter region). Results showed that, for rs321846600, a T>C mutation reduced the DNA–protein binding affinity (Fig 5E). For rs1111440035, both the M4-L-CC and M4-H-GG probes detected two common bands (Fig 7E, arrows 2 and 3 in lanes 2 and 5), whereas the M4-L-CC probe detected an additional shifted band (Fig 7E, arrow 1 in lane 5). When using U- M4-H-GG to compete with B- M4-L-CC, the shift was not eliminated (Fig 7E, arrow 1 in lane 8). Thus, the specific band suggested that unknown trans-acting factors could uniquely bind to the M4-L-CC probe and likely respond to the different promoter activities of BMP2 H and BMP2 L. However, for the M5-L-CC and M5-H-GG probes for rs80791204, no additional shifted bands were detected (Fig 7F), suggesting that these alleles are not different in binding affinity. Combined with the above results, this suggests that SNPs rs321846600 and rs1111440035 influence transcription factor binding and may alter the expression of the BMP2 gene, which thus impacts LMD trait expression. Discussion In this study, we pre-corrected LMD phenotypes and genotyped the DNA of 2,733 Yorkshire pigs using a 50k SNP chip for GWAS. Five SNPs in the genomic region between 15.51 and 16.31 Mb on SSC17 were significantly associated with LMD. To fully understand the GWAS results, we combined GWAS, LD block analysis, and Hi-C data to map the candidate genomic region related to LMD. Compared to traditional linkage analysis, LDLA is able to improve QTL detection and accurately map QTL locations [37,38]. Combined with the LDLA results, a ~10 kb QTL region (LD block2, 15.65 Mb to 15.75 Mb) and the BMP2 gene on chromosome 17 were associated with LMD in Yorkshire pigs. To further detect important functional variants in the transcriptional regulatory region, we used target region sequencing technology to increase the SNP density of the candidate genomic region. Together with SNP association analysis, obtaining information on cis-regulatory elements, and motif analysis, we were able to identify five variants were as important functional variant candidates. Dual-luciferase expression assays and EMSAs were used to evaluate these important candidate variants. Correspondingly, one SNP was detected that impacts enhancer activity and one SNP was detected that affects promoter activity of the BMP2 gene. Compared to traditional research methods which only used GWAS or integrated GWAS with gene expression data [39–41], our method provides conclusive evidence of a major QTL affecting LMD on chromosome 17 and further identifies candidate SNPs with functional relationships to LMD. Genome-wide association studies have been extensively used to identify genomic regions associated with important economic traits, however, revealing trait-related genetic mechanisms is still a research challenge. To identify the important candidate variants, 3D epigenomic characteristics were used to uncover important candidate functional sites as follows: 1) The TAD information surrounding the LMD-associated SNPs was used to dissect the SNP functions, which showed the genome region of significantly associated SNPs has limited opportunity to interact with genomic regions in other TADs. 2) The identified candidate genomic region (15.51 to 16.31 Mb) fully covers one sub-domain, and LD block 2 is located within this sub-domain. According to the property of sub-domains, the SNPs within LD block 2 have higher probability of interaction with genomic regions within the sub-domain. We admitted that we cannot totally exclude the possibility that variants outside block 2 may contribute to the LMD QTL effects, but the possibilities are much lower than the block 2. 3) The epigenetic data (enhancers and promoters) help to further identify candidate SNPs based on target re-sequencing data. 4) The dual-luciferase expression assays and EMSAs assist in validating the identified important variants. Overall, integrated analysis of genome-wide association studies and 3D epigenomic characteristics aid in understanding the function of variants governing LMD. Loin muscle depth is an important economic trait that influences the lean meat content of pigs. Revealing the genetic mechanism and identifying important markers contribute to the genetic improvement of LMD. LD analysis suggested that there were four blocks in the LMD significantly associated region, and these four blocks were found in the LMD-associated region and are related to body weight, body height, intramuscular fat content, average daily gain, backfat, and carcass length. LMD has been reported to be negatively correlated with backfat and positively correlated with carcass length [42,43]. Additionally, the most significant SNP in this LMD-associated region has been reported to be associated with LMA [44]. Therefore, the LMD-associated region may be the LMD QTL region. In our study, BMP2 was identified to be the candidate gene for LMD. BMP2, a member of the superfamily transforming growth factor beta (TGF-beta), has been reported to participate in adipogenesis [45–47], myogenesis [48,49], chondrogenesis, and osteogenesis [50–52]. BMPs, which are also known as activators or inhibitors, participate in specific stages of muscle progenitor cell differentiation during skeletal muscle development and regeneration [53–59]. Trait association analysis of human skeletal muscle volume revealed that BMP2 is associated with increased skeletal muscle volume [60]. In previous studies, BMP2 has been found to be related to carcass length [61,62], loin muscle area, body size, and several foot and leg (FL) structural soundness traits in pigs [44]. Endogenously or exogenously expressed BMP2 promotes adipogenesis [45,63]. Furthermore, BMP2 has been found to induce the upregulation of adipogenic gene expression, leading to increased intramuscular fat (IMF) deposition in castrated animals [64]. It has also been reported that LMD and LMA are negatively correlated with backfat depth and positively correlated with carcass length [42,43]. This evidence supports the hypothesis that BMP2 is a functional candidate gene regulating LMD in Yorkshire pigs. Overall, our work suggested that BMP2 expression plays a crucial role in driving LMD variation in Yorkshire pigs and that the C allele of rs321846600 reduces DNA–protein binding ability, which increased BMP2 enhancer activity. We identified two mutation sites in the promoter region of porcine BMP2, forming two haplotypes. The G alleles of rs1111440035 and rs80791204 are favorable as they can increase BMP2 promoter activity in cell reporter assays. Moreover, EMSA results indicated that the C allele at rs1111440035 binds to an unknown transcription factor, which leads to the C allele at rs1111440035 displaying weaker promoter activity than the G allele. To our knowledge, no transcriptional factors directly regulating BMP2 expression have previously been reported. In the future, we will perform supershift assays and protein pull-down experiments to characterize the transcription factors for SNP rs1111440035 (C>G) and rs321846600 (T>C) and explore the mechanism of regulation of BMP2 expression. Overall, our results indicate that both rs1111440035 and rs321846600 are likely important candidate mutations impacting LMD in Yorkshire pigs. Conclusion In our study, GWAS and epigenomics data, including ChIP-seq data and Hi-C data, were integrated to identify an LMD-related QTL region, as well as candidate genes within this region. A ~10 kb region (15659761 bp to 15755711 bp) on SSC17 was identified and determined to be the candidate LMD QTL region. The BMP2 gene was identified as the candidate gene most likely to be associated with LMD. SNPs rs1111440035 and rs321846600 are likely important candidate mutations impacting LMD in Yorkshire pigs. Our study made the first attempt to identify the major candidate genetic variants and candidate genes regulating an important production trait (LMD) in pigs through the integration of genome-wide association studies and 3D epigenomics. Supporting information S1 Table Description of gene information in the QTL region significantly associated with LMD. (XLS) Click here for additional data file. S2 Table Description of quantitative traits loci (QTL) in the regions significantly associated with LMD. (XLS) Click here for additional data file. S3 Table Summary of linkage disequilibrium and linkage analysis (LDLA). (XLS) Click here for additional data file. S4 Table Description of information for captured sequencing loci significantly associated with LMD. (XLS) Click here for additional data file. S1 Fig Motif analysis for the SNPs significantly associated with LMD by target regional sequencing. (TIF) Click here for additional data file. S1 Data The data of phenotype. (ZIP) Click here for additional data file. S2 Data The data of genotype (This file can be opened with Notepad++). (ZIP) Click here for additional data file. 10.1371/journal.pgen.1010820.r001 Author response to previous submission Submission Version0 10 Dec 2022 Attachment Submitted filename: reponse letter.docx Click here for additional data file. 10.1371/journal.pgen.1010820.r002 Decision Letter 0 Barsh Gregory S. Editor-in-Chief Groenen Martien Academic Editor © 2023 Barsh, Groenen 2023 Barsh, Groenen https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Submission Version0 11 Jan 2023 Dear Dr Xiang, Thank you very much for submitting your Research Article entitled 'Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs' to PLOS Genetics. The manuscript was fully evaluated at the editorial level and by independent peer reviewers. The reviewers appreciated the attention to an important problem, but raised some substantial concerns about the current manuscript. Based on the reviews, we will not be able to accept this version of the manuscript, but we would be willing to review a much-revised version. We cannot, of course, promise publication at that time. Should you decide to revise the manuscript for further consideration here, your revisions should address the specific points made by each reviewer. We will also require a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. If you decide to revise the manuscript for further consideration at PLOS Genetics, please aim to resubmit within the next 60 days, unless it will take extra time to address the concerns of the reviewers, in which case we would appreciate an expected resubmission date by email to plosgenetics@plos.org. If present, accompanying reviewer attachments are included with this email; please notify the journal office if any appear to be missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist. To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission. While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool.  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process. To resubmit, use the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder. We are sorry that we cannot be more positive about your manuscript at this stage. Please do not hesitate to contact us if you have any concerns or questions. Yours sincerely, Martien Groenen, PhD Academic Editor PLOS Genetics Gregory Barsh Editor-in-Chief PLOS Genetics While the manuscript provides clear evidence for a QTL affecting LMD, the results do not support the claim that two SNPs upstream of BMP2 are causative. Furthermore, I agree with reviewers 1 and 3 that the epigenetic data cannot be used to identify the causative variants within the region identified and therefore cannot be presented as a method for fine mapping QTLs. Both claims need to be toned down in the revised version of the manuscript Reviewer's Responses to Questions Comments to the Authors: Please note here if the review is uploaded as an attachment. Reviewer #1: This paper provides conclusive evidence for a major QTL affecting loin muscle depth (LMD) on chromosome 17 in a population of Yorkshire pigs and suggestive evidence for the functional significance for some SNPs located in putative regulatory regions of BMP2. The problem is that it is unknown if the candidate SNPs affect LMD or not. Therefore, the claim that two SNPs upstream of BMP2 are causative (L 49) is not justified. The authors do not have the genetic resolution to identify causal variants due to the broad QTL region (800 kb) with high LD among markers (Fig. 3). The initial mapping was done with a sparse SNP set, and the resequencing analysis carried out in this analysis show that the most strongly associated markers is located more than 200 kb away from the SNPs now claimed to be causative. The data are worth publishing but the authors need to tone down the arguments that these SNPs are causal it is more correct to state that these are candidate SNPs that may be functionally related to the LMD QTL. Also, I disagree with the statement that this paper presents a method for fine mapping QTLs. The epigenetic data used here are instructive for identifying functionally important regions of the genome but they cannot be used to exclude sequence variants from being causative and that is what is needed to fine map a QTL. This is the strength of genetic methods like linkage and LD mapping namely that they can exclude regions and thereby narrow down regions associated with a phenotype. What can be done as regards this locus would be to collect similar data from another pig population and also to search for sweep signals in this region. If this is a major locus under strong selection a clear sweep signal is expected among pigs selected for lean meat content. An alternative explanation is that this QTL is due to a haplotype effect caused by multiple closely linked QTLs within this interval so the effect of this region varies among pig populations. Specific comments Fig. 3. Please indicate coordinates for the region covered in Fig. 3A in order for readers to compare the regions mentioned in the text. Furthermore, clarify that all SNP data here are based on the sparse SNP chip. This figure illustrates the problem with the identification of causal variants because there is high LD across a large genomic region. The authors focus on the block2 region where BMP2 is located but SNPs in this region show complete or near complete LD with SNPs at the other end of the broad interval about 1 Mb away. Furthermore, it is confusing that the LRT peak is at around 9 Mb on chr 17 according to Fig. 3C whereas the text indicates that the peak is around 15.7 Mb?? L 298-305, resequencing analysis. The authors focus on the block 2 region 15.65-15.75 Mb, but this resequencing analysis (Table S4) shows that the most significantly associated SNPs are located in the interval 15.89-16.19Mb. The top SNP is located at 16.03 MB more than 280 kb from the candidate SNPs that are functionally evaluated. This is an important result that should be presented in a graph in one of the main figures. Since this is based on a dense SNP data set it provide information about the strength of association across the QTL region. L 320. I assume you should cite Fig. 5 here although Fig 5 shows results for only 5 regions and not 6 as indicated in the text. Furthermore, it is not the SNPs that can function as an enhancer, it is regions harboring the SNPs that may have enhancer function. Three constructs are used in Fig 5d: T, G and T-G. The difference between G and T-G is not clear to me and is not explained in the legend. Also, why is this not a candidate SNP when G vs T is significantly different? L 327-335. The description of the promoter polymorphism is a bit confusing. The authors refer to 5 SNPs but Fig. 6 suggests that the correct description is 4 SNPs + an InDel (5bp deletion or insertion), correct? Furthermore, the authors refer to rs numbers in the text but the figures use distance to TSS, these two annotations need to be connected in the text perhaps using the M1 to M5 annotations. Fig. 5e is confusing since constructs indicated below the bar is not annotated and do not seem to have been used so why include them in the figure?? Fig 7cd. There is much more comparative data available for mammals that can be used to assess conservation scores for these SNPs. The authors state that both M4 and M5 are evolutionary conserved but that is not correct for M5 because Fig. 7d shows that both C and G occur among other mammals. It would be better to refer these SNPs as REF and non-REF because these are polymorphism. For instance, for M4 the variant denoted MUT is the one present in all other mammals included here. Gel shift assays, Fig 5f and Fig. 7ef. Firstly, it is not possible to verify causative SNPs with gel shifts, it can merely support their functional significance but it cannot prove causality, change subtitle. The authors should add one more lane to these experiments and test if excess cold probe of the other allele can compete for binding or not. If a mutation inactivates binding, cold mutation probe will not be able to compete with the labeled oligo binding the nuclear factor. For instance, Fig 5f shows that Cold U-WT can compete with B-WT but it is of interest to know if cold U-MUT can compete with B-WT or not. Fig. 7e has a lot of background and is not convincing, it needs to be improved. The authors highlight three bands labeled 1, 2 and 3, they should explain their interpretation which of these are specific and which are not. L 367. This is not correct that you fine mapped the QTL you still have a large region, but you provided functional evaluation of candidate SNPs. L. 388-389 I disagree with this statement, it is important to distinguish genetic evidence and functional evidence (see main comment above). Minor comments L 68-70. It is not the number of SNPs that is most relevant here but the number of independent loci associated with the trait, because a single QTL may involve hundreds or thousands of SNPs that show statistical significance if whole genome sequencing has been performed. L 250. Us the same format consistently for genomic regions, so this should be indicated as 15.08 – 16.76 Mb L 254 I think it is better to write that “the LMD-related causal variants were expected…”. This is because you use genetic data to map causal variants but if these are regulatory they may act on distance and the genes regulated by the causal variants maybe located outside the interval. L. 265, change “mapping” to “map”, one example that some further improvement of English language is needed. Furthermore, should this be “QTLs” in pluralis? I assume you identified a single QTL that may include multiple causal variants but if so, these segregate as a single QTL. (However, the major QTL maybe composed of multiple closely linked QTLs but the authors do not have the genetic resolution to distinguish this architecture from a single QTL scenario). L 285 Delete “The” at the beginning of the sentence, I assume you have not yet characterized all cis-regulatory elements in the porcine genome. L 293 “two pig subpopulations” please clarify if these are selected from the same Yorkshire population used for the initial QTL mapping. L301-305 This sentence is incomplete and the language must be corrected. L 396, change to backfat L 420-421 It is more appropriate to write binds “an unknown transcription factor”. Table 2 is sorted by P-value but it is better to sort it based on genomic position since that gives a hint of the peak of association. L 686, one decimal place is sufficient here and replace “p-value” with “-log10(p)-value)” Fig. 2B. Why is this one single TAD-region, there is not much interaction between the 15.08 and 16.76 Mb Reviewer #2: Cograts to the authors for the good piece of scientific work showing the have identified not only a major gene, but also a QTN. In my opinion, all the cahnges suggested by the reviewers were accomplished and Engish is in good shape Reviewer #3: This paper describes a study in which a very strong statistical association was detected between a region on pig chromosome 17 and a key phenotype of commercial pigs, and subsequently, well-chosen genetic methods were used to dissect the functional differences between SNPs associated with the detected association. General comments Overall, this is a very interesting paper. The integration of functional analyses with the GWAS and LD-based analyses provide compelling results that implicate specific SNPs from the BMP2 gene in the association with loin muscle depth. My main concern, which could be addressed with rewording in various places, is that the title, abstract and conclusions state that 3D epigenomic characteristics helped them to dissect the basis of the genetic association but I do not see this. As I read the paper, the epigenomic characteristics did not really contribute to detecting the association or narrowing down the candidate SNPs, and rather, the GWAS and LD-based analyses were sufficient to detect the association. What I see as the novelty of the study is the effective use of the dual-luciferase and electrophoretic mobility shift assays to help identify which SNPs drive the strong association detected by GWAS/LDLA. Another issue in the manuscript is that the results based on target region sequencing are not well integrated with those for the SNP array results. As I understand it, the target region sequencing work was added in response to reviewers’ comments on the previous version of the manuscript and these results make an important contribution to the paper. However, the two sets of results need to be combined better to provide a more logical explanation of the study design and results. This applies to the Methods (161-179), Results (285-308) and Discussion (371-389) (some specific issues are also detailed below in Additional comments). There are some other sections where further details are required (see Additional comments, below). Finally, there are various grammar and wording problems. I have suggested some corrections, but further editing is required in some sections (see Grammar/wording, below). Additional comments Methods 107-135: I was confused by the treatment of phenotypes such that polygenic effects have been fitted twice (pre-correction + GWAS), using two different approaches. The authors should explain/justify this further. 153-160: The authors should provide further (brief) details here (the reader shouldn’t need to read the other paper to understand what has been done). 163-165: Seven out of how many SNPs? 164, 169, 170: Why 15.65-15.75 Mb on line 164, then subsequently 15.51-16.31 Mb? 167: Were the 732 individuals selected from high- and low-LMD groups or across the whole population? 167-169: Need to provide more details for the target region sequencing. 165, 172…: Which dataset was used for this association analysis described on 174-179? The 393 animals from high- and low-LMD or the 732 individuals? 180-202: I found this section confusing and further clarification is necessary. It was not clear to me which groups were being compared, i.e. What is meant by “experimental” and “control” groups (202)? Also, which pig(s) were used to generate the cells and provide DNA (189)? 207: As I understand, this was performed for 3 SNPs. Assuming so, this should be stated in the text (e.g. “For the three identified candidate SNPs…”) Results 244-246: More details are required regarding the TADs. The reader should not need to read Reference [18] in order to understand this paper. Also see wording comment below. 284-308: As mentioned above, these results are not integrated clearly with the results from the KASP genotyping. For example, why wasn’t there any overlap between the two sets of identified SNPs? More generally, the results of the two approaches should be linked as they are addressing the same aim. 353-362: This paragraph needs more details on how to interpret these results, in particular for rs1111440035 or rs80791204. For example, “no additional shifted bands were detected” (359-360), should be following by “suggesting that these alleles did not differ in binding affinity” (if that is correct?). I could not understand the results for rs1111440035 (357-360). Discussion 371-389: As mentioned above, this section needs editing to give a more logical explanation of the study design and interpretation of results. Grammar/wording 76: change “variations” to “variation” (2 places); change “these” to “this” 102: reword, e.g. “SNPs were mapped to the Sscrofa11.1 pig genome assembly.” 108: reword, the ssGBLUP analysis incorporates both genetic and non-genetic effects. 150: add “A” before “likelihood ratio test” 170-171: needs rewording (I don’t understand what is meant here by “meanwhile distributed”) 173: change to “R software” or “R” 180: change to “…promoter and enhancer regions” 182: change “predict” to “identify” 183: change to “affect a transcription factor binding site” 185-188: This section needs rewording --be more precise about their proximity of the two SNPs --I don’t follow why 5 fragments were cloned --be more precise about “Similar experiments…” --change “executed” to “carried out” 190: should this be “from genomic DNA”? 192: change “assay” to “assays” 194: change “cell” to “cells” 217-229: this section requires editing to correct grammar (e.g. verb tenses) and wording 234-236: change “standard errors” to “standard error” (3 places) 238: change to “A linear mixed model analysis was applied…” 245: need to reword “insulation of interactions” (I don’t have any idea what this means) 250: change to “in the same TAD region” 256: add “the” before “above-mentioned” 265-2666: change “performed the linkage disequilibrium (LD) analysis” to “investigated linkage disequilibrium (LD) patterns” 274: change “vital” to another word (“compelling”?) 287-288: change to “Based on whole genome sequencing data for 60 Yorkshire pigs” 288: remove “the” 289: change to “the pig genome” 290: should mention this is on SSC17 298-308: this section requires editing to correct various grammar and wording problems 320: change to “supporting the hypothesis that all 6 of these SNPs can function as enhancers” 322: change to “in the enhancer region” 323: change “rest” to “remaining” 340: change “affecting” to “affect” 344: change “differential” to “different” 357: add “the” before “T>C” 359: change “probe” to “probes” 360-362: change to “…this suggests that SNPs rs321846600 and rs1111440035 influence transcription factor binding and affect…” 392: remove “the” 392-393: change to “four blocks are found in the LMB-associated region” 396: typo (“backfat”) 397: do you mean “the most significant SNP”? 398-399: reword sentence; does “this region” refer to that reported in Reference 44? 413: change to “This evidence supports the hypothesis that …” 415-427: this section requires editing to correct grammar and improve wording ********** Have all data underlying the figures and results presented in the manuscript been provided? Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information. Reviewer #1: Yes Reviewer #2: Yes Reviewer #3: No: I presume that the journal will check that the data is made available if the paper is accepted. Currently the data is not available. ********** PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files. If you choose “no”, your identity will remain anonymous but your review may still be made public. Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy. Reviewer #1: No Reviewer #2: No Reviewer #3: No 10.1371/journal.pgen.1010820.r003 Author response to Decision Letter 0 Submission Version1 29 Mar 2023 Attachment Submitted filename: reponse letter.pdf Click here for additional data file. 10.1371/journal.pgen.1010820.r004 Decision Letter 1 Barsh Gregory S. Editor-in-Chief Groenen Martien Academic Editor © 2023 Barsh, Groenen 2023 Barsh, Groenen https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Submission Version1 7 May 2023 Dear Dr Xiang, Thank you very much for submitting your Research Article entitled 'Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs' to PLOS Genetics. The manuscript was fully evaluated at the editorial level and by independent peer reviewers. The reviewers appreciated the attention to an important topic but identified some concerns that we ask you address in a revised manuscript. We therefore ask you to modify the manuscript according to the review recommendations. Your revisions should address the specific points made by each reviewer. In addition we ask that you: 1) Provide a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. 2) Upload a Striking Image with a corresponding caption to accompany your manuscript if one is available (either a new image or an existing one from within your manuscript). If this image is judged to be suitable, it may be featured on our website. Images should ideally be high resolution, eye-catching, single panel square images. For examples, please browse our archive. If your image is from someone other than yourself, please ensure that the artist has read and agreed to the terms and conditions of the Creative Commons Attribution License. Note: we cannot publish copyrighted images. We hope to receive your revised manuscript within the next 30 days. If you anticipate any delay in its return, we would ask you to let us know the expected resubmission date by email to plosgenetics@plos.org. If present, accompanying reviewer attachments should be included with this email; please notify the journal office if any appear to be missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist. While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission. To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice. PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process. To resubmit, you will need to go to the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder. Please let us know if you have any questions while making these revisions. Yours sincerely, Martien Groenen, PhD Academic Editor PLOS Genetics Gregory Barsh Editor-in-Chief PLOS Genetics As reviewer 1 points out, you still oversell your data as providing conclusive evidence that the causal variants are located within a 10 kb region by using epigenetic data. Please tone down this statement and address the remaining comments made by the two reviewers. Reviewer's Responses to Questions Comments to the Authors: Please note here if the review is uploaded as an attachment. Reviewer #1: This paper reports convincing evidence for a major QTL on pig chromosome and report functional characterization of candidate SNPs in an enhancer and the promotor of the BMP2 gene. The data are good. However, the remaining problem is that the authors still oversell their data as providing conclusive evidence that the causal variants are located within a 10 kb region by using epigenetic data. The problem here is that the QTL region is large 800 Mb and functional data can neither prove or exclude the importance of a sequence variant for a genotype-phenotype relationship, this can only be achieved by genetic analysis. For instance, a functional assay such as a luciferase assay used here can show that a sequence variant affect function but it does not prove that this effect is important for the genotype-phenotype relationship. Further, functional assay may fail to reveal a significant effect of a sequence variant because the assay does not replicate the conditions when the sequence variant is important. For instance, in this study the authors use a kidney cell line to study the possible effect of sequence variants affecting muscle development (LMD). In conclusion, I think the paper is fine if the authors correct some remaining issues described below and are more realistic what they can and cannot conclude based on these data. Hopefully, my comments below are useful. Specific comments L302- (Narrowing down the QTL region). I am still not convinced that these data are sufficient to narrow down the QTL region from about 800 kb to about 10 kb. The reason is that the authors have used a sparse SNP panel. There could be other sequence variants within the 800 kb region with equally strong association to phenotype as those within LD block 2. For instance, Figure 3 shows that the SNP INRA0052824 that is located about 800 kb from the top SNP still shows D’=1 (or very close to 1) to the top SNP in the 10 kb region. Figure 4d. This figure illustrates very well that it is impossible to identify a 10 kb interval as the sole genomic region harboring causal mutations for the LMD QTL because the strong association with phenotype is over a broad region and may reflect a haplotype effect. To make this figure coherent I suggest that the authors mark their favorite region 15.65-15.75 in Fig. 4d. As far as I can see it is not at all the most strongly associated region in 4d. The authors should also explain in the legend that this plot tests for genetic association to the LMD phenotype. Figure 5. The text referring to these experiments is a bit confusing. The authors clone 800-bp genomic regions harboring 5 candidate SNPs. The problem is apparently that these 800 bp fragments may contain other SNPs in addition to the candidate SNP. This becomes apparent in Fig. 5d when the result first is highly significant but then the authors make the specific T to G change in the T construct and then there is no significant difference. The authors conclusion on line 364-365 “that SNP 15684170 is not an important site for enhancer activity” is wrong. This is still possible but the important conclusion is that the highly significant difference between the T and G construct must be caused by another SNP in this interval, perhaps that is a candidate causal sequence variant? The authors end this section by concluding that SNP rs321846600 is the best candidate SNP. However, it is not clear if the 800 bp fragment also carries other SNPs that may be important. The authors need to sequence the four 800 bp fragments and declare which sequence differences each region contains in a Supplementary Table. If the fragment with SNP rs321846600 contains multiple sequence variants they need to mutate the candidate SNP to prove that this is underlying the difference in luciferase activity between constructs. L432-445. I still find the authors arguments here problematic. They are mixing up genetic significance and functional significance. Epigenetic data cannot be used to narrow down a QTL region. It can be used to test the functional significance of SNPs within a QTL region. In this paper the authors have decided to focus on a 10 kb region but they have not excluded that sequence variants outside this 10 kb region but within the 800 kb region contribute to the QTL effect or even is more important than the possible effect of the variants within the 10 kb region. This paragraph needs to be modified to not be misleading to the field. In particular the sentence (The significant SNPs ..) on line 440-441 needs to be followed by something like this: “However, this does not exclude the possibility that sequence variants outside this region contribute to the LMD QTL effect in the genomic region 15.51-16.31 on SSC17.” L468-471 There is a conflict between these two studies. I recommend that “demonstrated” is replaced with “suggested” Minor comments L32: delete “the” L34: replace “the key” with “candidate” L45-46: change to “candidate SNPs that may be functionally related to” L51: change “major” to “candidate” L66. Change “loci” to “locus” L69. Change “the major” with “candidate” L277: Change “were” to “are” L346: Since these are predictions, I think this sentence is too strong and should read “Results suggested that these five SNPs may all disturb…” L357 and 359: delete “different carried” L396: change “supports the fact” with “shows” L397: change “conservative” to “conserved” Figure 5 and 6. Please add a few words to the legend and explain that PK15 is a porcine kidney cell line so that readers don’t have to search for this info in the M&M Figure 5e: It is unclear what WT and MUT refer to here. It would be much better to use the nucleotides instead alleles T and C. That would make it much easier to compare the results in Fig. 5c and 5e. It may also be better to change WT and MUT to H and L in Fig. 7 because it is not clear what is WT and MUT given that these are polymorphisms segregating in pig populations. Line 415: Add “using a 50k SNP chip” after “pigs” in order to explain that a sparse SNP panel has been used at this step. Reviewer #3: The authors have made substantial changes to the manuscript and addressed most of my concerns from the previous review. I think it will be suitable for publication following minor changes and editing, as described below. The authors state that they have had assistance from professional editors and/or native English speakers, however, the manuscript would benefit from further editing. There are many places where the language/grammar is problematic. Specific comments I am not convinced that including F_ST between American Yorkshire and American Landrace breeds is a good addition to the paper. In my view, this weakens what is otherwise a rigorous study. There is no reason to expect that these two breeds will be differentiated in this region of the genome, which has been identified using GWAS within a population of Yorkshire pigs. Or at least, the authors do not provide any justification for this analysis. In terms of the Results, they show a small peak overlapping the region identified by the GWAS but this may be coincidental. There is a higher peak in the 3’ direction (about which I also wouldn’t interpret too much). As they only looked at F_ST in a small part of the genome, there is no way to tell how extreme these peaks are on the genome-wide scale. I would recommend removing this section of the paper. Lines 278-293: I appreciate the addition of details from the previous study to which the authors refer. However, I found this section somewhat hard to follow. Are the authors saying that the Hi-C data helped them narrow down the region identified by GWAS (line 282)? I don’t see how this is the case if the TAD region is actually larger than the GWAS region. --line 279: change “isolate” to another word (“limit”?) --need to define “sequesters” --line 286 (and elsewhere in the paper): need to explain in the text what is meant by “interact” and “interactions” in this context --lines 291-293: please reword to clarify what is meant by “… SNPs are closely related to each other” and “… their target genes were expected in this TAD region” Figure 4d. The legend for this plot needs to be changed to explain that this shows results from an association analysis (as I understand it, as described on lines 199-206?). Line 348-365: Change “carried” to a different word (“cloned”?) Line 389: Change “NR2C1 cannot binding” (I don’t know what this means) Line 397: Change “conservative” to “conserved” Line 409: Change to either “do not differ” or “are not different” Line 422: Remove “newly developed” Line 426: Change “verify” to “evaluate” Line 430: Change “provides” to “identifies” Line 472: Change to “…increased BMP2 enhancer activity…” Line 478: Should this be “have previously been reported”? ********** Have all data underlying the figures and results presented in the manuscript been provided? Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information. Reviewer #1: Yes Reviewer #3: Yes ********** PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files. If you choose “no”, your identity will remain anonymous but your review may still be made public. Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy. Reviewer #1: No Reviewer #3: No 10.1371/journal.pgen.1010820.r005 Author response to Decision Letter 1 Submission Version2 5 Jun 2023 Attachment Submitted filename: response Letter.pdf Click here for additional data file. 10.1371/journal.pgen.1010820.r006 Decision Letter 2 Barsh Gregory S. Editor-in-Chief Groenen Martien Academic Editor © 2023 Barsh, Groenen 2023 Barsh, Groenen https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Submission Version2 7 Jun 2023 Dear Dr Xiang, We are pleased to inform you that your manuscript entitled "Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs" has been editorially accepted for publication in PLOS Genetics. Congratulations! Before your submission can be formally accepted and sent to production you will need to complete our formatting changes, which you will receive in a follow up email. Please be aware that it may take several days for you to receive this email; during this time no action is required by you. Please note: the accept date on your published article will reflect the date of this provisional acceptance, but your manuscript will not be scheduled for publication until the required changes have been made. Once your paper is formally accepted, an uncorrected proof of your manuscript will be published online ahead of the final version, unless you’ve already opted out via the online submission form. If, for any reason, you do not want an earlier version of your manuscript published online or are unsure if you have already indicated as such, please let the journal staff know immediately at plosgenetics@plos.org. In the meantime, please log into Editorial Manager at https://www.editorialmanager.com/pgenetics/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production and billing process. Note that PLOS requires an ORCID iD for all corresponding authors. Therefore, please ensure that you have an ORCID iD and that it is validated in Editorial Manager. To do this, go to ‘Update my Information’ (in the upper left-hand corner of the main menu), and click on the Fetch/Validate link next to the ORCID field.  This will take you to the ORCID site and allow you to create a new iD or authenticate a pre-existing iD in Editorial Manager. If you have a press-related query, or would like to know about making your underlying data available (as you will be aware, this is required for publication), please see the end of this email. If your institution or institutions have a press office, please notify them about your upcoming article at this point, to enable them to help maximise its impact. Inform journal staff as soon as possible if you are preparing a press release for your article and need a publication date. Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Genetics! Yours sincerely, Martien Groenen, PhD Academic Editor PLOS Genetics Gregory Barsh Editor-in-Chief PLOS Genetics www.plosgenetics.org Twitter: @PLOSGenetics ---------------------------------------------------- Comments from the reviewers (if applicable): ---------------------------------------------------- Data Deposition If you have submitted a Research Article or Front Matter that has associated data that are not suitable for deposition in a subject-specific public repository (such as GenBank or ArrayExpress), one way to make that data available is to deposit it in the Dryad Digital Repository. As you may recall, we ask all authors to agree to make data available; this is one way to achieve that. A full list of recommended repositories can be found on our website. The following link will take you to the Dryad record for your article, so you won't have to re‐enter its bibliographic information, and can upload your files directly:  http://datadryad.org/submit?journalID=pgenetics&manu=PGENETICS-D-22-01410R2 More information about depositing data in Dryad is available at http://www.datadryad.org/depositing. If you experience any difficulties in submitting your data, please contact help@datadryad.org for support. Additionally, please be aware that our data availability policy requires that all numerical data underlying display items are included with the submission, and you will need to provide this before we can formally accept your manuscript, if not already present. ---------------------------------------------------- Press Queries If you or your institution will be preparing press materials for this manuscript, or if you need to know your paper's publication date for media purposes, please inform the journal staff as soon as possible so that your submission can be scheduled accordingly. Your manuscript will remain under a strict press embargo until the publication date and time. This means an early version of your manuscript will not be published ahead of your final version. PLOS Genetics may also choose to issue a press release for your article. If there's anything the journal should know or you'd like more information, please get in touch via plosgenetics@plos.org. 10.1371/journal.pgen.1010820.r007 Acceptance letter Barsh Gregory S. Editor-in-Chief Groenen Martien Academic Editor © 2023 Barsh, Groenen 2023 Barsh, Groenen https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 15 Jun 2023 PGENETICS-D-22-01410R2 Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs Dear Dr Xiang, We are pleased to inform you that your manuscript entitled "Integrated analysis of genome-wide association studies and 3D epigenomic characteristics reveal the BMP2 gene regulating loin muscle depth in Yorkshire pigs" has been formally accepted for publication in PLOS Genetics! Your manuscript is now with our production department and you will be notified of the publication date in due course. The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript. Soon after your final files are uploaded, unless you have opted out or your manuscript is a front-matter piece, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers. Thank you again for supporting PLOS Genetics and open-access publishing. We are looking forward to publishing your work! With kind regards, Zsofia Freund PLOS Genetics On behalf of: The PLOS Genetics Team Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom plosgenetics@plos.org | +44 (0) 1223-442823 plosgenetics.org | Twitter: @PLOSGenetics ==== Refs References 1 Wang K , Liu D , Hernandez-Sanchez J , Chen J , Liu C , Wu Z , et al . Genome Wide Association Analysis Reveals New Production Trait Genes in a Male Duroc Population. PLoS One 2015, 10 (9 ):e0139207. doi: 10.1371/journal.pone.0139207 26418247 2 He Y , Ma J , Zhang F , Hou L , Chen H , Guo Y , et al . Multi-breed genome-wide association study reveals heterogeneous loci associated with loin eye area in pigs. J Appl Genet 2016, 57 (4 ):511–518. doi: 10.1007/s13353-016-0351-8 27183999 3 Wellcome Trust Case Control Consortium. Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature. 2007 Jun 7;447 (7145 ):661–78. doi: 10.1038/nature05911 17554300 4 Daetwyler HD , Capitan A , Pausch H , Stothard P , van Binsbergen R , Brondum RF , et al . Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. Nat Genet 2014, 46 (8 ):858–865. doi: 10.1038/ng.3034 25017103 5 Onteru SK , Gorbach DM , Young JM , Garrick DJ , Dekkers JC , Rothschild MF . Whole Genome Association Studies of Residual Feed Intake and Related Traits in the Pig. PloS one 2013, 8 (6 ):e61756. doi: 10.1371/journal.pone.0061756 23840294 6 Cho IC , Yoo CK , Lee JB , Jung EJ , Han SH , Lee SS , et al . Genome-wide QTL analysis of meat quality-related traits in a large F2 intercross between Landrace and Korean native pigs. Genet Sel Evol 2015, 47 (7 ):014–0080. 7 Debin NI , Zhao S , Xia X , Junyong HU , Liu W , Yunlong MA . Prediction of loin eye muscle area in live pig using loin eye muscle depth. Journal of Huazhong Agricultural University 2018. 8 Xue Y , Li C , Duan D , Wang M , Han X , Wang K , et al . Genome-wide association studies for growth-related traits in a crossbreed pig population. Anim Genet 2020. 9 Zhuang Z , Li S , Ding R , Yang M , Zheng E , Yang H , et al . Meta-analysis of genome-wide association studies for loin muscle area and loin muscle depth in two Duroc pig populations. PLoS One 2019, 14 (6 ):e0218263. doi: 10.1371/journal.pone.0218263 31188900 10 Thomsen H , Lee HK , Rothschild MF , Malek M , Dekkers JC . Characterization of quantitative trait loci for growth and meat quality in a cross between commercial breeds of swine. J Anim Sci 2004, 82 (8 ):2213–2228. doi: 10.2527/2004.8282213x 15318717 11 Grindflek E , Szyda J , Liu Z , Lien S . Detection of quantitative trait loci for meat quality in a commercial slaughter pig cross. Mamm Genome 2001, 12 (4 ):299–304. doi: 10.1007/s003350010278 11309662 12 Mohrmann M , Roehe R , Knap PW , Looft H , Plastow GS , Kalm E . Quantitative trait loci associated with AutoFOM grading characteristics, carcass cuts and chemical body composition during growth of Sus scrofa. Anim Genet 2006, 37 (5 ):435–443. doi: 10.1111/j.1365-2052.2006.01492.x 16978171 13 Tak YG , Farnham PJ . Making sense of GWAS: using epigenomics and genome engineering to understand the functional relevance of SNPs in non-coding regions of the human genome. Epigenetics Chromatin 2015, 8 (57 ):015–0050. doi: 10.1186/s13072-015-0050-4 26719772 14 Ernst J , Kheradpour P , Mikkelsen TS , Shoresh N , Ward LD , Epstein CB , et al . Mapping and analysis of chromatin state dynamics in nine human cell types. Nature 2011, 473 (7345 ):43–49. doi: 10.1038/nature09906 21441907 15 Maurano MT , Humbert R , Rynes E , Thurman RE , Haugen E , Wang H , et al . Systematic localization of common disease-associated variation in regulatory DNA. Science 2012, 337 (6099 ):1190–1195. doi: 10.1126/science.1222794 22955828 16 Claussnitzer M , Dankel SN , Kim KH , Quon G , Meuleman W , Haugen C , et al . FTO Obesity Variant Circuitry and Adipocyte Browning in Humans. The New England journal of medicine 2015, 373 (10 ):895–907. doi: 10.1056/NEJMoa1502214 26287746 17 Miao Y , Mei Q , Fu C , Liao M , Liu Y , Xu X , et al . Genome-wide association and transcriptome studies identify candidate genes and pathways for feed conversion ratio in pigs. BMC genomics 2021, 22 (1 ):294. doi: 10.1186/s12864-021-07570-w 33888058 18 Zhao Y , Hou Y , Xu Y , Luan Y , Zhou H , Qi X , et al . A compendium and comparative epigenomics analysis of cis-regulatory elements in the pig genome. Nature communications 2021, 12 (1 ):2217. doi: 10.1038/s41467-021-22448-x 33850120 19 Legarra A , Aguilar I , Misztal I . A relationship matrix including full pedigree and genomic information. J Dairy Sci 2009, 92 (9 ):4656–4663. doi: 10.3168/jds.2009-2061 19700729 20 Christensen OF , Lund MS . Genomic prediction when some animals are not genotyped. Genet Sel Evol 2010, 42 :2 . 21 Madsen P , Jensen J . A user’s guide to DMU. A package for analyzing multivariate mixed models. Version 6, release 5.2. University of Aarhus. Center for Quantitative Genetics and Genomics Dep of Molecular Biology and Genetics, Research Centre Foulum, Tjele, Denmark 2013. 22 Yin L , Zhang H , Tang Z , Xu J , Yin D , Zhang Z , et al . rMVP: A Memory-efficient, Visualization-enhanced, and Parallel-accelerated tool for Genome-Wide Association Study. bioRxiv 2020:2020.2008.2020.258491. doi: 10.1016/j.gpb.2020.10.007 33662620 23 VanRaden PM . Efficient methods to compute genomic predictions. Journal of dairy science 2008, 91 (11 ):4414–4423. doi: 10.3168/jds.2007-0980 18946147 24 Barrett JC , Fry B , Maller J , Daly MJ . Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 2005, 21 (2 ):263–265. doi: 10.1093/bioinformatics/bth457 15297300 25 Zhang F , Zhang Z , Yan X , Chen H , Zhang W , Hong Y , et al . Genome-wide association studies for hematological traits in Chinese Sutai pigs. BMC Genet 2014, 15 :41. doi: 10.1186/1471-2156-15-41 24674592 26 Legarra A , Fernando RL . Linear models for joint association and linkage QTL mapping. Genet Sel Evol 2009, 41 :43. doi: 10.1186/1297-9686-41-43 19788745 27 Grindflek E , Lien S , Hamland H , Hansen MH , Kent M , van Son M , et al . Large scale genome-wide association and LDLA mapping study identifies QTLs for boar taint and related sex steroids. BMC Genomics 2011, 12 :362. doi: 10.1186/1471-2164-12-362 21752240 28 Druet T , Georges M . A hidden markov model combining linkage and linkage disequilibrium information for haplotype reconstruction and quantitative trait locus fine mapping. Genetics 2010, 184 (3 ):789–798. doi: 10.1534/genetics.109.108431 20008575 29 Sartelet A , Druet T , Michaux C , Fasquelle C , Geron S , Tamma N , et al . A splice site variant in the bovine RNF11 gene compromises growth and regulation of the inflammatory response. PLoS Genet 2012, 8 (3 ):e1002581. doi: 10.1371/journal.pgen.1002581 22438830 30 Garcia-Gamez E , Gutierrez-Gil B , Sanchez JP , Arranz JJ . Replication and refinement of a quantitative trait locus influencing milk protein percentage on ovine chromosome 3. Anim Genet 2012, 43 (5 ):636–641. doi: 10.1111/j.1365-2052.2011.02294.x 22497507 31 Lander ES , Botstein D . Mapping mendelian factors underlying quantitative traits using RFLP linkage maps. Genetics 1989, 121 (1 ):185–199. doi: 10.1093/genetics/121.1.185 2563713 32 Qiao R , Gao J , Zhang Z , Li L , Xie X , Fan Y , et al . Genome-wide association analyses reveal significant loci and strong candidate genes for growth and fatness traits in two pig populations. Genet Sel Evol 2015, 47 :17. doi: 10.1186/s12711-015-0089-5 25885760 33 Durand NC , Robinson JT , Shamim MS , Machol I , Mesirov JP , Lander ES , et al . Juicebox Provides a Visualization System for Hi-C Contact Maps with Unlimited Zoom. Cell systems 2016, 3 (1 ):99–101. doi: 10.1016/j.cels.2015.07.012 27467250 34 Quinlan AR , Hall IM . BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 2010, 26 (6 ):841–842. doi: 10.1093/bioinformatics/btq033 20110278 35 Robinson JT , Thorvaldsdóttir H , Turner D , Mesirov JP . igv.js: an embeddable JavaScript implementation of the Integrative Genomics Viewer (IGV). bioRxiv 2020:2020.2005.2003.075499. doi: 10.1093/bioinformatics/btac830 36562559 36 Castro-Mondragon JA , Riudavets-Puig R , Rauluseviciute I , Berhanu Lemma R , Turchi L , Blanc-Mathieu R , et al . A JASPAR 2022: the 9th release of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 2022 Jan 7;50 (D1 ):D165–D173.34850907 37 Farnir F , Grisart B , Coppieters W , Riquet J , Berzi P , Cambisano Net al. Simultaneous mining of linkage and linkage disequilibrium to fine map quantitative trait loci in outbred half-sib pedigrees: revisiting the location of a quantitative trait locus with major effect on milk production on bovine chromosome 14. Genetics 2002, 161 (1 ):275–287.12019241 38 Meuwissen TH , Karlsen A , Lien S , Olsaker I , Goddard ME . Fine mapping of a quantitative trait locus for twinning rate using combined linkage and linkage disequilibrium mapping. Genetics 2002, 161 (1 ):373–379. doi: 10.1093/genetics/161.1.373 12019251 39 Godia M , Reverter A , Gonzalez-Prendes R , Ramayo-Caldas Y , Castello A , Rodriguez-Gil JE , et al . A systems biology framework integrating GWAS and RNA-seq to shed light on the molecular basis of sperm quality in swine. Genet Sel Evol 2020, 52 (1 ):72. doi: 10.1186/s12711-020-00592-0 33292187 40 Luo W , Xu J , Li Z , Xu H , Lin S , Wang J , et al . Genome-Wide Association Study and Transcriptome Analysis Provide New Insights into the White/Red Earlobe Color Formation in Chicken. Cellular physiology and biochemistry: international journal of experimental cellular physiology, biochemistry, and pharmacology 2018 , 46 (5 ):1768–1778. doi: 10.1159/000489361 29705805 41 Yan Z , Huang H , Freebern E , Santos DJA , Dai D , Si J , et al . Integrating RNA-Seq with GWAS reveals novel insights into the molecular mechanism underpinning ketosis in cattle. BMC genomics 2020, 21 (1 ):020–06909. doi: 10.1186/s12864-020-06909-z 32680461 42 Wilson KB , Overholt MF , Hogan EK , Schwab C , Shull CM , Ellis M , et al . Predicting pork loin chop yield using carcass and loin characteristics. J Anim Sci 2016, 94 (11 ):4903–4910. doi: 10.2527/jas.2016-0610 27898928 43 Holland LA , Hazel LN . RELATIONSHIP OF LIVE MEASUREMENTS AND CARCASS CHARACTERISTICS OF SWINE. Journal of Animal Science 1958(3 ):825. 44 Fan B , Onteru SK , Du ZQ , Garrick DJ , Stalder KJ , Rothschild MF . Genome-wide association study identifies Loci for body composition and structural soundness traits in pigs. PloS one 2011, 6 (2 ):e14726. doi: 10.1371/journal.pone.0014726 21383979 45 Guiu-Jurado E , Unthan M , Böhler N , Kern M , Landgraf K , Dietrich A , et al . Bone morphogenetic protein 2 (BMP2) may contribute to partition of energy storage into visceral and subcutaneous fat depots. Obesity 2016, 24 (10 ):2092–2100. doi: 10.1002/oby.21571 27515773 46 Ji X , Chen D , Xu C , Harris SE , Mundy GR , Yoneda T . Patterns of gene expression associated with BMP-2-induced osteoblast and adipocyte differentiation of mesenchymal progenitor cell 3T3-F442A. Journal of bone and mineral metabolism 2000, 18 (3 ):132–139. doi: 10.1007/s007740050103 10783846 47 Sottile V , Seuwen K . Bone morphogenetic protein-2 stimulates adipogenic differentiation of mesenchymal precursor cells in synergy with BRL 49653 (rosiglitazone). FEBS letters 2000, 475 (3 ):201–204. doi: 10.1016/s0014-5793(00)01655-0 10869556 48 Katagiri T , Akiyama S , Namiki M , Komaki M , Yamaguchi A , Rosen V , et al . Bone morphogenetic protein-2 inhibits terminal differentiation of myogenic cells by suppressing the transcriptional activity of MyoD and myogenin. Experimental cell research 1997, 230 (2 ):342–351. doi: 10.1006/excr.1996.3432 9024793 49 Katagiri Takenobu , Yamaguchi Akira . Bone morphogenetic protein-2 converts the differentiation pathway of C2C12 myoblasts into the. Journal of Cell Biology 1994 , 127 (6 ):1755–1755. 50 Salazar VS , Gamer LW , Rosen V . BMP signalling in skeletal development, disease and repair. Nature reviews Endocrinology 2016, 12 (4 ):203–221. doi: 10.1038/nrendo.2016.12 26893264 51 Zhou N , Li Q , Lin X , Hu N , Liao JY , Lin LB , et al . BMP2 induces chondrogenic differentiation, osteogenic differentiation and endochondral ossification in stem cells. Cell and tissue research 2016, 366 (1 ):101–111. doi: 10.1007/s00441-016-2403-0 27083447 52 Lv Y , Gao CW , Liu B , Wang HY , Wang HP . BMP-2 combined with salvianolic acid B promotes cardiomyocyte differentiation of rat bone marrow mesenchymal stem cells. The Kaohsiung journal of medical sciences 2017, 33 (10 ):477–485. doi: 10.1016/j.kjms.2017.06.006 28962818 53 Cossu G , Borello U . Wnt signaling and the activation of myogenesis in mammals. The EMBO journal 1999, 18 (24 ):6867–6872. doi: 10.1093/emboj/18.24.6867 10601008 54 Dietrich S , Schubert FR , Healy C , Sharpe PT , Lumsden A . Specification of the hypaxial musculature. Development 1998, 125 (12 ):2235–2249. doi: 10.1242/dev.125.12.2235 9584123 55 Linker C , Lesbros C , Stark MR , Marcelle C . Intrinsic signals regulate the initial steps of myogenesis in vertebrates. Development 2003, 130 (20 ):4797–4807. doi: 10.1242/dev.00688 12917295 56 Tajbakhsh S. Stem cells to tissue: molecular, cellular and anatomical heterogeneity in skeletal muscle. Current opinion in genetics & development 2003, 13 (4 ):413–422. doi: 10.1016/s0959-437x(03)00090-x 12888016 57 Fazzalari NL . Bone fracture and bone fracture repair. Osteoporosis international: a journal established as result of cooperation between the European Foundation for Osteoporosis and the National Osteoporosis Foundation of the USA 2011 , 22 (6 ):2003–2006. doi: 10.1007/s00198-011-1611-4 21523400 58 Sartori R , Sandri M . BMPs and the muscle-bone connection. Bone 2015, 80 :37–42. doi: 10.1016/j.bone.2015.05.023 26036170 59 Sartori R , Schirwis E , Blaauw B , Bortolanza S , Zhao J , Enzo E , et al . BMP signaling controls muscle mass. Nature genetics 2013, 45 (11 ):1309–1318. doi: 10.1038/ng.2772 24076600 60 Devaney JM , Tosi LL , Fritz DT , Gordish-Dressman HA , Jiang S , Orkunoglu-Suer FE , et al . Differences in fat and muscle mass associated with a functional human polymorphism in a post-transcriptional BMP2 gene regulatory element. Journal of cellular biochemistry 2009, 107 (6 ):1073–1082. doi: 10.1002/jcb.22209 19492344 61 Blaj I , Tetens J , Preuss S , Bennewitz J , Thaller G . Genome-wide association studies and meta-analysis uncovers new candidate genes for growth and carcass traits in pigs. PloS one 2018, 13 (10 ):e0205576. doi: 10.1371/journal.pone.0205576 30308042 62 Falker-Gieske C , Blaj I , Preuß S , Bennewitz J , Thaller G , Tetens J . GWAS for Meat and Carcass Traits Using Imputed Sequence Level Genotypes in Pooled F2-Designs in Pigs. G3 (Bethesda). 2019 Sep 4;9 (9 ):2823–2834. doi: 10.1534/g3.119.400452 31296617 63 Denton NF , Eghleilib M , Al-Sharifi S , Todorčević M , Neville MJ , Loh N , et al . Bone morphogenetic protein 2 is a depot-specific regulator of human adipogenesis. Int J Obes 2019, 43 (12 ):2458–2468. 64 Jung DJS , Baik M . Up-regulation of bone morphogenetic protein and its signaling molecules following castration of bulls and their association with intramuscular fat content in Korean cattle. Scientific reports 2019, 9 (1 ):19807. doi: 10.1038/s41598-019-56439-2 31875043