Detection of CNVs in wrinkled skin cases and controlsBased on the Illumina HiSeq2500 platform, we performed whole-genome resequencing and genome-wide CNVs detection on 20 cases of wrinkled skin (WSC) and 63 controls from the Xiang pig population. The mapped read depth ranged from 9.76 to 18.45, with an average depth of 13.75 per sample (Supplementary Table S1). The CNVs were called by CNVnator and CNVcaller software from the NGS data based on the Read Depth method. A total of 7893 CNVRs were identified in all samples, including 5753 losses and 2140 gains (Supplementary Tables S2-3). The CNVRs covered 47.34% of the porcine genome (Sscrofa 11.1) in all populations. The group of wrinkled skin cases had 5562 CNVRs, with 4559 losses and 1003 gains. These CNVs were located on all chromosomes except for Y chromosome with a mean size of 147,229.16 bp ranging from 600 to 1,996,600 bp (Table 1; Fig. 1). All the CNVR maps for 83 pigs were showed in Fig. 2. In addition, sequence coverage did not correlate with the number of CNVs identified in each individual (Fig. 3).
Table 1 Descriptive statistics of copy number variants identified for two groups.Fig. 1Size distribution of CNVRs. The distribution of CNV fragment size in the groups with wrinkled skin (cases) and those without (controls) is presented herewith. A consistent trend is evident in the distribution of CNV size between the two groups, with small fragments of CNV being more prevalent in the loss type and large fragments of CNV being more prevalent in the gain type.Fig. 2A comprehensive chromosome mapping of all CNVRs. Two types of CNVR were identified: gain (green) and loss (grey). The appearance of vertical bars on a chromosome indicates the position of discrete types of CNVs.Fig. 3The relationship between the number of variants and the sequence coverage for each sample. No correlation was observed between sequence coverage and the number of CNVRs identified in each individual.Comparison with CNVRs identified in previous reportsComparing our CNV data to those reported previously in several other studies for pigs, and the results showed that our data has a 27.26% overlap in CNV events with other researchers (Table 2). Comparison showed that the average overlap rate was is 9.91%, which was is higher than the average overlap rate of 2.13% found in previous studies18. The comparison also showed that the overlapping CNV events accounted for 61.93% of the total events, while the newly identified CNVRs accounted for 38.07%. These results suggest that the detected CNV events were highly reliable. Here we adopted the principle that two CNVRs share at least one base to identify overlapping CNV events.
Table 2 Comparison between CNVR detected in this study and CNVR reported previously.Population stratification analysisTo discover CNVR that were differentiated between WSC and controls, the pairwise F-statistic (Fst) between wrinkled skin and cases controls samples were calculated by using the NGS data at each locus. As shown in Supplementary Table S3 (Fig. 4), pairwise Fst values showed generally large (Fst = 0.151–0.248) to very large (Fst = 0.25–0.44) genetic divergence between WSC and controls. About 86.87% (n = 6819) of the CNVRs showed nonsignificant genetic differentiation (Fst = 0.0–0.05), while 11.96% (n = 944) of the CNVRs showed moderate genetic differentiation (Fst = 0.05–0.15), 0.988% (n = 77) and 0.203% (n = 16) of the CNVRs showed large (Fst = 0.15–0.248) to very large (Fst = 0.25–0.44) genetic differentiation, respectively. A Fst value of 0.15 was used as a threshold for identifying potentially selective CNVR. In total, we identified 93 WSC-controls stratified CNVRs, including 17 gains and 76 losses, which had a length ≥ 0.8 kb, and the largest event was CNVR_4865, a 1.715 Mb gain in 11 chromosome. In order to study the potential function of the population-stratified CNVRs, the CNV coverage area were annotated. As shown in Table S4, there were 87 known genes (such as VCAN, YKT6, CAMK2B, GCK, WNT9A, FOXO1, THSD1), 57 novel genes, 1 pseudogene, and 45 noncoding RNA genes overlapped with 62 CNVRs. A total of 31 of the 93 CNVRs were located in the intergenic region, without covering any genes.Fig. 4Identification of regions with differentiation in WSC compared to controls. A total of 93 differential CNVRs were identified using a threshold of Fst > = 0.15 (indicated by the red line).We further investigated the impact of population-stratified CNVRs on covering genes by Variant Effect Predictor (VEP). CNVR_1233 deletion (loss) indicated a high impact predicted to be loss-of-function, including feature truncation, transcript ablation, transcript amplification and frame shift variants on host genes. Total of 7 DUPs (gains) were characterized to be whole-gene duplication which generally caused the copy gain of an entire gene. CNVR_4865, for example, overlaped 11 ensembl genes, including SOX11, FOXO1, MRPS31, SLC25A15, THSD1, CKAP2, NEK3, NEK5, ALG11, CCDC70 and one novel gene. Whole-gene duplication results in transcript amplification, which affects gene expression in a dose-dependent manner.GO and KEGG analysis of the genes covered by population-stratified CNVRsIn order to gain a deeper understanding of the biological functions of the genes overlapped with the population-stratified CNVRs and the mechanism of the relationship between CNVRs and disease, we conducted gene ontology (GO) and KEGG enrichment analysis on the 87 known genes using the KOBAS tool. There were 47 genes that mapped to 149 KEGG pathways (Supplementary Table S5), including 54 significant pathways (P < 0.05), which mainly enriched in Insulin and Insulin resistance, FoxO, Notch, N-Glycan biosynthesis, Amino sugar and nucleotide sugar metabolism, Mitophagy – animal, Cancer, Th1 and Th2 cell differentiation, Endocrine resistance, AGE-RAGE signaling pathway in diabetic complications, and ErbB signaling pathway. There were 85 genes mainly enriched in 63 significant GO terms (P < 0.01), including protein N-linked glycosylation, extracellular matrix, metalloendopeptidase activity, mitophagy, elastic fiber assembly, negative regulation of apoptotic process, acetylglucosaminyltransferase activity, glutathione peroxidase activity, response to mitochondrial depolarization, extracellular matrix organization, hyaluronic acid binding and others.Screening of candidate CNVRs related with wrinkled skin casesWe then investigated the CNVR-covering genes that involved in the biological processes and pathways of skin aging, such as apoptosis, survival, endocytosis, cell differentiation, oxidative stress, DNA damage and repair, extracellular matrix, immune and metabolism, and combined with the previously reported genes and databases (https://www.ncbi.nlm.nih.gov/gene/?term=aging) related to aging, we identified 22 CNVRs putatively associated with wrinkled skin cases, containing 61 known genes (Table 3). Three CNVRs (CNVR_1195, CNVR_1854, CNVR_4865) were predicted to cause transcript amplification effects on genes (such as FOXO1, ADAMTS2, MAPK9, SQSTM1, NEK5, NEK3). Ten known genes, including VCAN, CAMK2B, TIMP1, and LGALS2, were affected by 9 deletion variants resulting in transcript ablation, stop lost, feature truncation. Ten other CNVRs may impact the intronic regions of the AGK, SDR1, DMD, LPGAT1, and ADAMTS18 genes.
Table 3 Descriptions of the significant CNVRs associated with wrinkled skin cases.CNV assaySix RNA-seq datasets were used to investigate the impact of candidate CNVRs on the expression of overlapping genes in case and control skin samples. Of the 61 known genes overlapping 22 candidate CNVRs, 56 were expressed in the skin (Supplementary Table S6). Wrinkled skin cases and controls showed differential expression of twelve genes (Fig. 5), including TIMP1, FOXO1 and ARAF. The expression levels of TIMP1, FOXO1, and ARAF genes were 1.7–3 times lower in the case group compared to the control group.Fig. 5The expression trends of 12 CNVR-overlapping differential genes in cases of wrinkled skin (WSC) and control skin (XP).To further clarify the functional significance of the candidate CNVRs in skin wrinkling traits in aromatic pigs, we validated CNV_VCAN (CNVR_1233), CNV_TIMP1 (CNVR_7451). We evaluated the copy numbers of CNV_VCAN, CNV_TIMP1 in F2 populations of the WSC × LW family using quantitative qPCR. All individuals were divided into three categories according to 2−ΔΔCT values (Supplementary Table S7): those with 2−ΔΔCT values ranging from 1.5 to 2.5 as normal type and those with values less 1.5 or more than 2.5 represented loss and gain, separately. The types of normal, loss and gain were detected in the 58 samples from the F2 population. The frequency of copy number variation of the three CNV genes in the F2 population was inconsistent. The proportion of losses accounted for by CNV_TIMP1 was 39.66%, followed by CNV_VCAN at 58.62%. CNVs may affect the phenotype by altering the transcription of genes through a dosage effect. We evaluated the relationship between CNV types of CNV_VCAN, CNV_TIMP1 and skin wrinkles traits in the F2 population using a general linear model (Table 4). The skin wrinkles trait was expressed as wrinkled fraction, which is equal to the wrinkle value divided by the profile value. The study suggests that skin wrinkling was affected by different CNV genotypes, specifically CNV_TIMP1. Increasing the copy number of CNV_TIMP1 resulted in more wrinkles on the F2 population (P < 0.05). RNA-seq data also showed a significant increase in the expression of gained CNV_TIMP1 copy number phenotypes at 6 and 12 months of age, particularly at 12 months (P < 0.05). However, the correlation between CNV_VCAN copy number and skin wrinkles was not strong. However, at 12 months of age, the skin wrinkle index was significantly higher in the loss type compared to the normal type.
Table 4 Association analysis of CNV types with wrinkled skin traits.
Genome-wide identification of copy number variations in wrinkled skin cases of Xiang pigs
