Integration of multi-omics approaches for functional characterization of muscle related selective sweep genes in Nanchukmacdon

Pig as a food source serves daily dietary demand to a wide population around the world. Preference of meat depends on various factors with muscle play the central role. In this regards, selective breeding abled us to develop “Nanchukmacdon” a pig breeds with an enhanced variety of meat and high fertility rate. To identify genomic regions under selection we performed whole-genome resequencing, transcriptome, and whole-genome bisulfite sequencing from Nanchukmacdon muscles samples and used published data for three other breeds such as Landrace, Duroc, Jeju native pig and analyzed the functional characterization of candidate genes. In this study, we present a comprehensive approach to identify candidate genes by using multi-omics approaches. We performed two different methods XP-EHH, XP-CLR to identify traces of artificial selection for traits of economic importance. Moreover, RNAseq analysis was done to identify differentially expressed genes in the crossed breed population. Several genes (UGT8, ZGRF1, NDUFA10, EBF3, ELN, UBE2L6, NCALD, MELK, SERP2, GDPD5, and FHL2) were identified as selective sweep and differentially expressed in muscles related pathways. Furthermore, nucleotide diversity analysis revealed low genetic diversity in Nanchukmacdon for identified genes in comparison to related breeds and whole-genome bisulfite sequencing data shows the critical role of DNA methylation pattern in identified genes that leads to enhanced variety of meat. This work demonstrates a way to identify the molecular signature and lays a foundation for future genomic enabled pig breeding.

Pig is the most studied animal model to date with Mitochondrial DNA (mtDNA) analysis tracked down its origin from Eurasian wild-boar 1 . Among other animals, pig share a deep connection with human civilization and played a critical role in fulfilling the feed demands 2 . Estimates for the consumption of animal food suggest that pork demand will rise to 36% of the overall meat consumption by 2025 and hence require attention to sustain and improve meat quality with its production to secure global food demand 3 .
The Korean peninsula is one of the largest pig consuming country with its Jeju native black pig (JNP) is an indigenous variety of Korean pig with high-quality meat content, redness and nutrition value 4,5 . Superior taste leads to an increase in JNP demand with every passing day but less productivity and fertility made it difficult to sustain the demand 6,7 . Due to huge demand but less productivity intrusion of different breeds leads to diversion of industry focus on alternative economic viable options 8,9 . These imported pig breeds possessed the excellent genetic potential for high production and their growth have been reported to more than 0.5 kg of weight per day but limited with meat quality content 10,11 . This ultimately threatened the indigenous variety of pig with almost reached extinction until the government put serious effort and involved in saving the native pig breed by close monitoring the growth and use.
To address the issue and sustain the demand of JNP in the entire region with food security for the long run, the emphasis has been given to developing a breed with high productivity and meat quality. An in-house breeding program was started to develop a breed with indigenous pig features and have a high fertility rate. Marker-based multiple inter-crosses using a strict selection of breeding pigs (Jeju Native black pig, Duroc, Landrace) accelerated the generation of outstanding progeny containing high meat-quality breed. In the course of continuously close monitoring and breeding program using modern biological methods, Nanchukmacdon a mixed breed was developed which maintained superior characteristics features in generations. A blind test for sensory reflexes to Nanchukmacdon was performed for meat quality and it was voted similar or better values than those for bacon meats. The carcasses of these newly developed black pigs (Nanchukmacdon) showed the significantly higher levels of intramuscular fat deposition (p < 0.05), the redness and the yellowness (p < 0.05) but not in the lightness (p > 0.05) in meat color contents 12 .
Identification of genomic regions undergone positive selection is a potent approach to delineate genes that help in adaptation to environmental factors and responsible for the phenotypic diversity. In the last decade, many GWAS studies have been conducted in this regards and supported by statistical advancement analysis to pin down significant results from the driven data. These approaches already helped in the identification of different genomic regions with selection signals, suggesting the contribution of the region in influencing certain characteristics related to phenotypic or genotypic composition in different breeds. To further extent, genes identified from whole-genome sequencing (WGS) may have evolved to adapt to the conditions but identifying the one expressing and controlling the trait features ideally govern these characteristics and identifying such differentially expressed genes for different trait could help us in identifying markers of interest to develop breeds with such high potential. In this search, RNAseq and whole-genome bisulfite sequencing (WGBS) analysis approaches are the attractive approach for the identification of differentially expressed genes and analysis of methylation role and have been used in various studies regarding such analysis [13][14][15] .
In this study, we have evaluated the genetic closeness of Nanchukmacdon with the related species and identified the selective sweep genes that allow the enhancement in the characteristics features possessed by the breed. We presented an unbiased approach combining WGS from different breeds of pig, RNAseq data from muscles of closely related species amongst them. We used statistically established methods such as cross-population extended haplotype homozygosity (XP-EHH) 16 and cross-population composite likelihood ratio (XP-CLR) 17 statistics to detect selection signatures from closely related breeds and Nucleotide Diversity analysis was performed in Nanchukmacdon using vcftools 18 . Finally, WGBS data analyzed for Nanchukmacdon to observe the methylation pattern in identified genes.

Results
Population structure analyses. PCA plot analysis describes the separation of species in 2 dimensional view. From the present analysis on 43 resequencing and RNA-seq sample analysis, we observed the clear association of 4 species. Nanchukmacdon is a breed between duroc, Landrace, and JNP and from the plot, it is conclusive that the genetic constituent is closer to duroc and with a certain contribution of JNP (Fig. 1a,b).
Positive selective signature in Nanchukmacdon population. In total, we have obtained 1154, 1296, and 1666 putative selection regions with p-values < 0.05 in XP-EHH and top 1% score limited to 574, 745, and 675 in XP-CLR putative positively selected genes test statistics in Nanchukmacdon from the three breeds Duroc, JNP, and Landrace respectively (Additional_data 1). Based on the analysis, we further narrow down the obtained results by overlapping XP-EHH and XP-CLR results and limited to 37, 41, and 39 genes in the statistical analysis.
Identification and analysis of differentially expressed genes (DEGs) in muscle tissue. DESeq was used to identify statistically significant differences in gene expression obtained by featurecounts 18 . A cutoff value of fold change ≥ 1 and adjusted FDR correction p-value < 0.05 was selected to obtain DEGs between different breeds NC_JNP_Muscles, NC_DU_Muscles, and NC_LR_Muscles (Additional_data 2). Amongst identified DEGs 1655 were found common in all the breeds (Fig. 2a). The overall relationship between different breed was depicted by Volcano Plot (Fig. 2b).
Gene ontology and functional profiling studies. Separate analysis was performed for identified positive selective sweep genes identified from XP-EHH and XP-CLR score and DEGs in Nanchukmacdon from Duroc, JNP, and landrace. The functional annotations of genes were categorized into three groups such as molecular function, cellular component, and biological process. Most significant (corrected p-value < 0.05) GO terms such as regulation of synaptic plasticity (   www.nature.com/scientificreports/ Fig S1a). KEGG pathway analysis reveals the involvement of major pathways varies from Metabolic pathway, Fatty acid biosynthesis, ErbB signalling pathway, Adipocytokine signalling pathway, Calcium signalling pathway and Oxidative phosphorylation are some (Fig. 2c, Supplementary Fig. S2). Manhattan plot for muscles WGS data analysis was generated for NC with DU, JNP, and LR by XP-EHH and XP-CLR score and annotated with commonly identified differentially expressed genes with CMplot 19 . XP-EHH and XP-CLR score was plotted against the genomic position with the autosomal chromosomes in different colors (Fig. 3a).
Validation of identified genes. All the identified genes were incorporated in the string database for protein-protein interaction analysis and it was observed that they share different modules of networks in governing important characteristics features (Fig. 4a) 20 . The muscle tissue from Nanchukmacdon was collected and used for expression analysis for seven randomly selected genes. The transcripts selected for validation are SERP2, UGT8, MELK, ZGRF1, FHL2, NCALD, UBE2L6 as mentioned in Additional_table 3. The RT-PCR experiment for selected samples was executed with three replicates for each sample. GAPDH and Beta-actin were used as an endogenous control for normalizing the quantification cycle (Cq) value. The contrast is depicted in (Fig. 4b). Results indicate that the expression of selected genes to control is significantly correlated with high confidence p-value = 0.00232 and R = 0.931. WGBS data analysis was performed for Nanchukmacdon to analyze the methylation in identified genes. As expected from the observed pattern, DNA methylation level sharply decreased near 5 kb upstream region of TSSs and dropped to the lowest outside the region (Fig. 4c), methylation level remains stable after promoter region contributing to structural stability and regulation of gene expression. CpG Island was less expressed inside than outside of 5 kb CpG Island (Fig. 4d) 21 . Individual methylation pattern for all the identified genes confirms the pattern of no methylation corresponding with the distribution of gene promoters, usually prone to transcription ( Supplementary Fig. S3).  www.nature.com/scientificreports/

Discussion
Genomic selection has been the main tool supported by statistical analysis in genetic improvement of economically important traits 22 . Selective signature ideally helps in stabilizing traits that make a breed unique with its features. Studying these traits allows a better understanding of breed and help in developing enhanced breed with high economic values. In this regards, Nanchukmacdon, a mixed breed was developed that exhibits exceptional meat quality with a high fertility rate and performed selective sweep analysis to understand the genes that triggered the meat characteristics. WGS data has been used for various GWAS related studies [23][24][25][26] , Incorporation of RNAseq analysis approach to WGS data provide us with a more detailed and better understanding to the identified genes and how they behave in the present system which makes specific trait different than the parental generations. Recently, various studies have been published reporting selective genes in different breeds over a period of time with main focus on identifying selective signature associated with the breeds and how these genes expressing in the system [27][28][29] . As meat quality is accessed in terms of the carcass, feed conversion efficacy, color, taste, juiciness and majorly consist of 75% of the muscles tissue 30 . In this regards, to identify selective sweep genes in closely related breeds we used statistically established methods such as cross-population extended haplotype homozygosity (XP-EHH) 16 and cross population composite likelihood ratio (XP-CLR) 17 statistics in order to detect selection signatures from closely related breeds; two approaches were used as each has its own advantages. XP-EHH compares haplotype lengths of populations to detect selective sweeps when the allele has approached or achieved fixation in one population but remains polymorphic in the other population. XP-CLR is a statistic based on allele frequency differentiation across populations providing an advantage to detect older signals and selection on standing variation. Subsequently, we integrated RNAseq analysis data of muscles to identify the muscles associated genes in the mixed breed with their closed associated pig varieties. www.nature.com/scientificreports/ After identifying common genes amongst different varieties exhibiting positive selective signature identified using XP-EHH and XPCLR statistical test, identification of common genes expressing in muscles tissue of Nanchukmacdon limited the total number of genes amongst different breeds to 11 genes with 4 each in JNP, DU, and LR respectively. Amongst UGT8 was identified in DU and LR. Here, JNP 4 (MELK, SERP2, GDPD5, FHL2) genes, DU (UGT8, ZGRF1, NDUFA10, EBF3) genes and in LR (UGT8, ELN, UBE2L6, NCALD) genes were identified [Table 1]. Similarly, identification of negatively expressed genes in Nanchukmacdon has been performed and common genes amongst JNP, DU, and LR with the cutoff parameter for log2fold change > 1 were identified to be limited with 13 genes, Additional_table 1.
Our findings on the genetic relationship amongst different breeds agree with the previously reported studies on different pig breeds using WGS data 31 , signifying the grouping of Nanchukmacdon with Duroc, Jeju native pig, and Landrace (Fig. 1a-c) with visualization of DEGs were performed using heatmap analysis (Fig. 1d). The individuals from the Nanchukmacdon, landrace, duroc and JNP population were grouped according to their origin as identified by PCA (Fig. 1a,b). Our results indicate that Nanchukmacdon is phylogenetically closer to landrace (Fig. 1c). To estimate individual ancestry, admixture proportions were assessed without defined population information using ADMIXTURE with The individual population was grouped into separate clusters at K = 4 with the lowest cross validation error ( Fig. 1e and Supplementary Fig. S4) and these results were also confirmed from unrooted tree (Fig. 1c). In selection signature analysis, we further performed a comparative analysis with closely related breeds. Our study reveals a series of well-known and novel genes reported in muscles related to biological processes and their high expression pattern after getting positively selected in the mixed breed. RNAseq methodology allowed us to further understand the changes held in selective sweep genes during the evolution. We have seen major expression alteration in mineral content or related pathways by the different expression pattern of selective sweep genes specifically muscles related (FHL2, EBF3), calcium ion channel route (NCALD, UBE2l6), metabolism, and fatty acid metabolism pathways (UGT8, GDPD5, MELK) when it compares to their closely related breeds. The genes under selection were further investigated using Nucleotide diversity approach for the selective sweep genes. Low nucleotide diversity signifies the stabilization and activation of genes under specific circumstances 32,33 and in Nanchukmacdon with their respective paired breed and identified the low diversity reported over the genomic coordinates and haplotype diversity in comparison to all other breeds (Fig. 3b, Supplementary Figs. S5 and S6).
Amongst identified genes, UDP Glycosyltransferase 8 (UGT8) was found to be highly expressed in Nanchukmacdon w.r.t. Duroc and Landrace variety of pig. Although there is limited understanding of UGT8 role with meat quality previous studies relates its involvement in either lipid metabolism, sphingolipid metabolism,  40 . Predicted functional partner from string database also confirmed a close association with epidermal growth factor receptor that could affect the morphology of the muscles 41 . Expression of muscles related genes significantly influenced by the expression pattern of FHL2 and found to control various genes such as MyoD1, MyH3 and MyoG that play a central role in the development of muscles. Similarly, NCALD positively selected gene from NC_LR presents in intracellular cellular component and observed for their role in calcium mediated signaling. It is a calcium sensor which directly interacts regulates actin and clathrin. Similarly, UBE2L6 was reported in molecular functions related to ATP-binding, energy metabolism and nucleotide-binding. It is involved in ubiquitination of multiple substrates 42,43 and direct involvement in obese related pathways 44 . GDPD5 or glycerophosphodiester phosphodiesterase 2 (GDE2) was identified as positively selective sweep genes whose functional annotation and string database information for protein-protein interaction analysis resulted in the limited lead for direct involvement in any pathway or function but suggested a close association with ACSL3 gene which play a central role in fatty acid oxidation 45 . ELN identified as the positively expressed selective gene found in the extracellular matrix that provides structural support, biochemical or biomechanical cues for cells or tissues by structure lying external to one or more cells 46,47 . ZGRF1 is not well characterized and function is not known in Sscrofa but the gene is reported to be associated with translation, transcription, nonsense-mediated mRNA decay, RNA decay, miRNA processing, RISC assembly, and pre-mRNA splicing. EBF3: Positive selection of EBF3 leads to an up-regulation of myogenic regulatory factors including MyoD and Myf5 required for muscles development 48 . The function of Ebf3 outside of the neuronal system, however, has limited understanding. EBF3 play role in muscles specific transcription and have critical role in relaxation by directly regulating the expression of a Ca 2+ pump 49 . SERP2 (stress associated endoplasmic reticulum protein family member 2), a positively selected differentially expressed gene with limited understanding of any direct role in muscles related pathway or any process in pig is reported to be differentially expressed in splay leg piglets 50 .

Conclusion
It is a well-established fact that artificial selection has greatly shaped pig genomic variability during the process of domestication. The variation developed during the event helps the local industry to proliferate and fulfil the local meat demands efficiently. Our primary objective of this study was to identify selective sweep genes that were expressed in muscles and we were able to limit 13 potent genes with their important role in the muscles building process. GO analysis showed various pathways vary from regulation of synaptic plasticity, clathrindependent endocytosis, and positive regulation of neuron differentiation were significantly enriched. Similarly, KEGG pathway showed that Metabolic pathways, Calcium signaling pathway, and Endocytosis were significantly enriched in selective sweep genes. These results provide a better understanding of the role of identified genes in regulating muscles related pathways and further information to genomic evolution and selective mechanism which could help develop an enhanced breed with high muscles content.

Methods
To identify selective sweep genes in Nanchukmacdon we have performed multiple analysis and a flowchart has been developed for better understanding of the work (Fig. 5).
Sampling and data collection. This study aimed to find selection signatures of selective sweep genes using WGS data and integrating RNAseq analysis approaches in Nanchukmacdon, and other pig breeds (JNP, Landrace, and Duroc) to identify differentially expressed selective sweep genes and their role in the biological process. Samples were taken from healthy Male pigs belong to the same farm in Jeju Island with an average age of 2 years for (N = 1-6), and 3 years for (N = 7-10) (Additional_Table 4). All the experimental procedures were verified and approved by the National Institute of Animal Science, and carried out in compliance with the ARRIVE guidelines 51 . Whole-genome re-sequencing was performed from the blood sample taken from postharvest Nanchukmacdon (N = 10), RNAseq data was generated for Nanchukmacdon (N = 5) pair-end data after isolation of muscle tissue using TRIzol method following the manufacturer guideline and reported earlier 15 . Similarly, gDNA from Nanchukmacdon muscles was subjected to bisulfite conversion using the fragment size (250 bp ± 25 bp), WGBS was performed with MethylMiner Methylated DNA Enrichment kit, and then a sequencing library was constructed using the Illumina Paired-end sequencing on an Illumina, NovaSeq, 150bpX2. Whole genome bisulfite sequencing (N = 5) was performed from the taken sample. To filter variants and avoid possible false positives, the "VariantFiltration" argument of GATK was adopted with the following options: SNPs with MQ (mapping quality) > 40.0, MQRankSum < − 12.5, ReadPosRankSum < − 8.0 www.nature.com/scientificreports/ and quality depth (unfiltered depth of non-reference samples; low scores are indicative of false positives and artifacts) < 2.0 were filtered 59 . BEAGLE version 4.1 60 was used to infer the haplotype phase and impute missing alleles for the entire set of swine populations simultaneously. After all the filtering processes, a total of ~ 26 million SNPs were retained and used for further analysis.
Phylogenetic tree, admixture and principal component analysis. To accurately describe the population of the crossbred pig, we used SNP data from different breeds and employed SNPRelate R package to perform principal component analysis 61 . Subsequently, Newick file was prepared and viewed in ggtree 62,63 .
Results were further analyzed for the distribution of breeds in different coordinates. Similarly, the input file was implemented using unsupervised based clustering method by a program ADMIXTURE to estimate the breed composition of individual animals 64 . The analysis was run with K (number of breeds '4') ranging from 2 to 4 to depict the genetic background of Nanchukmacdon and graphical display of the output was performed in R.
Detection of genomic regions with putative signals of selection. Using the whole SNP sets defined from NC, DU, LR and JNP, the method cross-population extended haplotype homozygosity (XP-EHH) 16 and cross-population composite likelihood ratio method (XP-CLR) 17 was used to detect genome-wide selective sweep regions. XP-EHH assesses haplotype differences between two populations and is designed to detect alleles that have increased in frequency to the point of fixation or near fixation in one of the two populations being compared 65 . Whereas, XP-CLR is based on the linked allele frequency difference between two populations and is a unidirectional method to find the pattern with regional allelic frequency difference in-between population 66 . Developments in RNA-seq technology enable more comprehensive investigation of the transcriptome for gene expression studies 67 . The statistical analysis is also critical in transcriptomic studies using RNA-seq; specifically, for the normalization of quantitative measurements of expression 68,69 and detection of DEGs 70 .
The PE reads were checked for the quality assessment using FastQC 56 and removed low quality reads by Trimmomatic 71 using parameters leading:3 trailing:3 slidingwindow:4:15 headcrop:13 minlen:36 before proceeding sequence alignment. All quality-filtered PE reads were aligned to Sscrofa genome (Sscrofa11.1) from the University of California Santa Cruz (UCSC) 72 using Hisat2 73 and reads were counted using FeatureCount 18 . DESeq2 was used to identify differentially expressed genes 74 .
Gene ontology analysis detection and annotation of candidate genes. Lists of differentially expressed genes with FDR < 0.05 in Nanchukmacdon w.r.t. Duroc, JNP, and Landrace were compiled and submitted to DAVID v6.8 server 36 for functional annotation and enrichment analysis. For each list, enriched Gene Ontology (GO) 75 Biological Processes, Molecular functions and Cellular Compartments. These terms were then clustered semantically using the ReviGO server 76 . Enriched functions throughout the whole transcriptome of Nanchukmacdon with elevated GO-term function and the clustered lower-level GO-terms. The letter corresponds to letters found in the treemap for Biological process, Molecular function and Cellular compartment (Additional_data 3, Additional_table 2, Supplementary Fig. S1b). The functional annotation of commonly identified selective sweep genes with differentially expressed genes was performed with the Database for Annotation, Visualization and Integrated Discovery (DAVID) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The genomic coordinates of the regions with high XP-EHH and XP-LCR score for 10k window with 10k bin size were used as input data and used R package BiomaRt 77,78 to fetch the coordinates against sscrofa11.1 database for getting the gene_id information of the respective regions with highly significant results 78 . Subsequently, the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (https:// david. ncifc rf. gov/ 36,79 was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway 80 and Gene Ontology (GO) 75 enrichment analyses. The GO terms and KEGG pathways with corrected p-value < 0.05 were considered significant. REVIGO 76 and Clusterprofiler R package 81 were used for summarizing the GO terms.
WGBS data analysis. The analysis for WGBS data was performed using reproducible genomics analysis pipeline PiGx-bsseq to understand methylation patterns in identified genes 82 . Where sequence was initially performed for a quality check using trim_galore 83 and alignment were subjected to the filtration of duplicate reads with sam_blaster and sorted using SAMtools 54 afterwards mapped to the reference genome of sscrofa11.1 using Bismark 84 . Bismark methyl extractor was performed to measure the methylation in CpG context. Validation. WGBS was performed for Nanchukmacdon muscles tissue to see methylation pattern in the identified genes and gene wise methylation visualization was performed to observe methylation level at CpG island by importing the CpG methylation file extracted from bismark methylation extract into SeqMonk visualization tool. The coordinates for identified genes were fetched from Ensembl and each gene was visualized for methylation and found correlating results (Supplementary Fig. S3).
Trizol method (Invitrogen, UK) was used for the total RNA isolation. Qubit fluorometer (Invitrogen, UK), NanoDrop (Thermo Scientific, USA) and Bioanalyzer (Agilent, UK) were used for analyzing the quality of the isolated total RNA. High Capacity cDNA Reverse Transcription Kit (Applied Biosystems™, 4368814) was used for synthesizing the cDNA and RT-PCR was performed by SYBR Green Realtime PCR Master Mix (TOYOBO, QPK-201T).
Ethics approval and consent to participate. In this study, N refers to number of animals and All the experimental procedures were verified and approved by the Ethics committee of National Institute of Animal with ethical approval no: NIAS20181295. www.nature.com/scientificreports/