Comprehensive Genome and Transcriptome Analysis Identifies SLCO3A1 Associated with Aggressive Behavior in Pigs

Abstract

Copy number variation (CNV) represents a significant reservoir of genetic diversity within the genome and exhibits a strong association with economically valuable traits in livestock. The manifestation of aggressive behavior in pigs has detrimental effects on production efficiency, immune competency, and meat quality. Nevertheless, the impact of CNV on the aggressive behavior of pigs remains elusive. In this investigation, we employed an integrated analysis of genome and transcriptome data to investigate the interplay between CNV, gene expression changes, and indicators of aggressive behavior in weaned pigs. Specifically, a subset of pigs comprising the most aggressive pigs (MAP, n = 12) and the least aggressive pigs (LAP, n = 11) was purposefully selected from a herd of 500 weaned pigs following a mixing procedure based on their composite aggressive score (CAS). Subsequently, we thoroughly analyzed copy number variation regions (CNVRs) across the entire genome using next-generation sequencing techniques, ultimately revealing the presence of 6869 CNVRs. Using genome-wide association study (GWAS) analysis and evaluating variance-stabilizing transformation (V_ST) values, we successfully identified distinct CNVRs that distinguished the MAP and LAP counterparts. Among the prioritized CNVRs, CNVR-4962 (designated as the top-ranked p-value and V_ST value, No. 1) was located within the Solute Carrier Organic Anion Transporter Family Member 3A1 (SLCO3A1) gene. The results of our analyses indicated a significantly higher (p < 0.05) copy number of SLCO3A1 in the MAP compared to the LAP. Furthermore, this increased copy number exhibited a positive correlation with the CAS of the pigs (p < 0.05). Furthermore, we integrated genomic data with transcriptomic data from the temporal lobe to facilitate the examination of expression quantitative trait loci (eQTL). Importantly, these observations were consistent with the mRNA expression pattern of SLCO3A1 in the temporal lobe of both MAP and LAP (p < 0.05). Consequently, our findings strongly suggest that CNVs affecting SLCO3A1 may influence gene expression through a dosage effect. These results highlight the potential of SLCO3A1 as a candidate gene associated with aggressive traits in pig breeding programs.

Keywords: CNV; SLCO3A1; aggressive behavior; animal welfare; dosage effect.

Figure 1

The characteristics and descriptive statistics of CNVRs. (A) The distribution of CNVRs and variation types of the MAP and LAP in autosomes and the X chromosome, where “Low” to “High” shows the gene density on the pig chromosomes. The yellow square, green circle and purple triangle represent the both type, gain type and loss type, respectively. (B) The length distribution of the CNVRs. The red, green, blue and purple colors represent all type, gain type, loss type and both type, respectively. (C) The distribution of positions and frequency on the genome of CNVRs. (D) The CNVR quantity was positively correlated with the length of the chromosomes (p < 0.05). (E) The total length of the fragment of CNVRs on each chromosome was proportional to the CNVRs quantity (p < 0.01). (F) The length of CNVRs in the chromosomes was not significantly statistically related to the length of the chromosomes (p > 0.05). (G) The sequencing depth of sequencing data was not significantly correlated with the CNVRs quantity (p > 0.05).

Figure 2

Performed enrichment analysis, GWAS analysis and V_ST analysis of CNVRs. (A,B) are the GO enrichment analysis and KEGG pathway analysis of all genes annotated in CNVRs, respectively. Bubbles of different colors represent different p-values. (C) GWAS analysis was used to identify CNVRs with significant effects on pig aggression. Among them, CNVR-4962 was ranked first with a p-value of 4.42 $\times$ 10⁻⁷. (D) Evaluate the differential CNVRs based on the V_ST value between the MAP and the LAP. Among them, including the CNVR-4962 (No. 1).

Figure 3

QTLs overlapped analysis of the CNVRs with the top 1% V_ST value.

Figure 4

Analysis of transcriptomic data in the temporal lobes. (A) PCA plot of the TL-MAP and TL-LAP. TL-LAP and TL-MAP represent the temporal lobe of the LAP and the temporal lobe of the MAP, respectively. (B) The heat map shows the overall picture of gene expression in the TL-MAP and TL-LAP. (C) The volcano plot shows that there are 366 up-regulated genes and 452 down-regulated genes in the MAP relative to the LAP. (D,E) GSEA analysis showed that genes upregulated in the MAP were associated with the “Glycerolipid metabolism” pathway, while genes downregulated in the MAP were associated with the “Tryptophan metabolism” pathway. (F–H) GSVA analysis showed that the “Steroid biosynthesis” gene set was increased in the MAP, and the gene set related to “Tryptophan metabolism” upward in the LAP, TL-LAP and TL-MAP represent temporal lobe of the LAP and temporal lobe of the MAP, respectively.

Figure 5

Cluster analysis of DEGs in the TL−MAP and TL−LAP. (A) Keans method divides the DEGs in the TL−MAP and TL−LAP into seven clusters. (B) The Mfuzz method divides the DEGs in the TL-MAP and TL-LAP into nine clusters. Images (C,D) represent GO enrichment analysis of genes highly expressed in the TL-MAP (clusters 5, 6 and 7 in the Kmeans method, and clusters 1, 4 and 7 in Mfuzz mothed), where “BP”, “CC” and “MF” represent the biological process, cellular component and molecular function, respectively. Images (E,F) represent KEGG pathway enrichment analysis of genes highly expressed in the TL-MAP.

Figure 6

Integrated analysis of the genomic and transcriptomic data. (A) 407 trans-eQTL genes are associated with the “Alzheimer disease” and “Neuroactive ligand-receptor interaction”. (B) Interaction network of 3 CNVRs and 5 neighboring genes with eQTL effects; the yellow line indicates positive regulation, and the purple line indicates negative regulation. Graphs (C,D) represent the relative expression in the SLCO3A1 in the TL-MAP and TL-LAP, and (C,D) are based on transcriptomic data and qPCR mothed, respectively. Images (E,F) are the evaluation results of the copy number of SLCO3A1 via CNVcaller and qPCR (via qPCR: n = 228). Images (G,H) are the relationship between the copy number of SLCO3A1 and CAS. (I) ROC curves of SLCO3A1 copy number as a genetic marker in the MAP and LAP. (J) Single-cell transcriptome data showed extensive expression of SLCO3A1 in the frontal, parietal, temporal and parietal lobes, and in the hypothalamus. * represents p < 0.05, ** represents p < 0.01.