Summary of the phenotypes and genotypes
The basic statistics of 42-days body weight, subcutaneous fat weight, subcutaneous fat percentage, abdominal fat weight, and abdominal fat percentage for each population were shown in Additional file 2. The heritability, genetic correlations, and phenotypic correlations of each trait were shown in Table 1. The estimates of heritability range from 0.42 to 0.61. Trait SFW presented the highest heritability (0.61) estimates, while the AFP was the lowest (0.42). The highest phenotypic correlation (0.95) and genetic correlation (0.98) were shown between AFW and AFP, and the lowest phenotypic correlation (0.33) and genetic correlation (0.26) were calculated between BW42 and AFP. The phenotypic and genetic correlations among these fat-related traits (SFW, SFP, AFW, and AFP) were all greater than or equal to 0.62.
Table 1
|
BW42
|
SFW
|
SFP
|
AFW
|
AFP
|
BW42 (g)
|
0.43
|
0.79
|
0.44
|
0.57
|
0.33
|
SFW (g)
|
0.87
|
0.61
|
0.84
|
0.80
|
0.62
|
SFP (%)
|
0.50
|
0.94
|
0.49
|
0.76
|
0.71
|
AFW (g)
|
0.67
|
0.88
|
0.87
|
0.49
|
0.95
|
AFP (%)
|
0.26
|
0.70
|
0.84
|
0.98
|
0.42
|
*Heritability (bold, diagonal), genetic correlations (below diagonal), and phenotypic correlations (above diagonal) for body weight and fat traits |
Abbreviations: BW42, body weight at day 42; SFW: skin fat weight; SFP: skin fat percentage; AFW: abdominal fat weight; AFP: abdominal fat percentage |
The population of Pop2014, Pop2019, and Pop2020 were obtained 1778Gb, 4799.83Gb, and 4709.25 Gb clean data, respectively. The average reads mapping alignment rate was 96.57%. The variant calling pipeline identified 14,220,859 SNPs within the duck genomes. After quality filtration and exclusion of variants of uncharacterized chromosomes, 8,448,069 SNPs and 573,457 independent SNPs passed our filters. A principal component analysis (PCA) was constructed to examine the relatedness among three populations (Additional file 3).
Genome-Wide Association Study identified the candidate variants and genes
A total of 648 chromosome-wide significant SNPs (1/573457, P < 1.74E − 06) and 242 genome-wide significant SNPs (0.05/573457, P < 8.72E − 08) across 40 chromosomes were identified by loci-based analysis. A list of all the genome-level significantly associated variants is shown in Additional file 4 and a list of all the chromosome-level significantly associated variants is shown in Additional file 5. There were 71, 47, 155, 64, and 14 genome-wide significant SNPs identified for traits BW42, SFW, SFP, AFW, and AFP, respectively. In total, 39 candidate genes were annotated from 242 genome-wide significant SNPs within 5Kb. The estimated genomic inflation factor λ ranged from 1.027 to 1.066 among 5 traits (Additional file 6), suggesting no population stratification in the studied population after specifying Sex, Birthday, and Batch effects as covariates. Manhattan and Q-Q plot of association results from genome-wide association analysis of each trait are shown in Additional file 6. We identified a 7Mb region between nucleotide position 54 Mb and 61 Mb in chromosome 4 that showed a significant association with the 42-days body weight phenotype with the most significant of 3.86E-10 (rs458603786, T > C, intergenic variant, P = 3.86E-10). As for fat traits, two genome-wide significant SNPs (rs117362452, G > A, synonymous variant, and rs291299411, G > C, intron variant) were detected to be correlated with SFW, SFP, AFW, and AFP, which were located at 7.36 Mb on Chr11 and 1.29 Mb on Chr 29. There was a high genetic correlation between subcutaneous fat and abdominal fat, 11 genes (Additional file 5) were common genes between subcutaneous and abdominal fat were identified among a total of 35 significant genes from SFW, SFP, AFW, and AFP. By querying the functional annotations of these genes, we found that CLN8 is a potential gene affecting lipid synthesis and transport, which was related to SFW, SFP, and AFW. The joint Manhattan plot of SFW and SFP showed in Fig. 1A.
ATAC-seq analysis of mutations in the upstream of CLN8
The five SNPs (rs322493594, rs322493619, rs322493641, rs322493648, rs322493651) were located in 1920 bp upstream of the CLN8 gene on chromosome 3. The rs322493651 was the most significant SNP accounting for 1.46% of the genetic variance, and 11.180% of the phenotypic variance. Linkage analysis revealed that they were in high linkage disequilibrium (LD, r2 > 0.9; Fig. 1B). These results prompted us to speculate that the causative mutation(s) should be the five variations. To confirm the association of the five mutations with the subcutaneous fat phenotype, 616 ducks were genotyped for the CLN8 heterozygous mutation and 3 for the homozygous, the wild type, and heterozygous mutation of haplotype with 5 SNPs significantly decreased the subcutaneous fat weight (P < 0.0001) (Fig. 1C). We used the ATAC-seq to investigate the chromatin accessibility landscape of subcutaneous preadipocytes. Fragments from the nucleosome-free regions (NFR) are expected to be enriched around the transcription start site (TSS) of genes (Additional file 7). Transcription factor binding sites were identified from 22,493,347 bp to 22,494,147 bp, which is located in the upstream of gene CLN8 (Fig. 2A). To exploit motif information, the JASPAR database was used to predict the motif binding site[37]. Different motif binding site was predicted as shown in Fig. 2B between wild-type sequence (up) and mutation sequence (down). ATAC-seq analysis showed that transcription factor binding sites were identified in a region close to the haplotype.
Active region screening of the CLN8 promoter
Based on Pekin duck genomic DNA, we designed primers in the proposed promoter region of CLN8 and amplified seven deletion fragments of different lengths: 3009 bp, 2519 bp, 2026 bp, 1349 bp, 811 bp, 501 bp, and 285 bp (Fig. 3A). To explore the active region of the CLN8 promoter, seven truncated fragments of the CLN8 promoter were amplified and inserted into the pGL3.1-Basic vector (Promega, USA) using restriction enzymes Kpn I and Xho I (Fig. 3B), and were co-transfected into ICP1 cells with the pRL-TK vector for luciferase activity detection. The results showed that the luciferase activity of pGL3.1-F3, pGL3.1-F2, and pGL3.1-F1 were significantly increased compared to that of pGL3.1-Basic, which was the control group (P < 0.05), the promoter activity of pGL3.1-F4 was significantly reduced (P < 0.05), and the activity of pGL3.1-F3 was significantly increased compared to that of pGL3.1-F4 (P < 0.05) (Fig. 3C). These results demonstrated that the region from − 1884 to -1207 bp upstream of the transcription start site of CLN8 was the promoter core transcription active region.
The functional polymorphism SNPs decreased the transcriptional activity of CLN8
We next investigated the effect of the five SNPs (Wild-type: -A-T-C-C-G- and Mutant-type: -T-C-A-A-C-) variant on the transcriptional activity of the CLN8 gene. Wild and Mutant haplotype CLN8 plasmids were verified before transfection by direct sequencing (Fig. 4A). To assess the effects of the polymorphisms, we generated two luciferase reporter gene constructs that share identical backbone sequences except for the polymorphisms, as shown in Fig. 4B, reporter gene expression driven by the A-T-C-C-G-containing CLN8 promoter (Wild-type) was greater than that driven by the T-C-A-A-C-containing counterpart (Mutant-type) (P < 0.001). This result demonstrated that significantly lower transcriptional activity of the mutant haplotype was observed when compared with the wild types.
CLN8 promotes the differentiation of avian adipocytes
To further determine the roles of CLN8 in avian adipocyte differentiation, we first performed gain-of-function experiments by using piggyback delivery of CLN8 into ICPs, and we detected its expression in the overexpression group and control group cells. As shown in Figs. 5A and 5B, a transfected CLN8 clone can significantly increase CLN8 expression at the level of transcription and translation by qPCR and western blot. Further, we detected the expression of adipocyte markers, PPARγ, and FABP4, and the results showed that they were significantly increased in CLN8OE cells (Fig. 5C, 5D). As shown by oil red O staining of neutral lipids and detection of the relative lipid content, over-expression of CLN8 in these cells (CLN8OE) significantly facilitated lipid formation (Fig. 5E, 5F). These results indicate that CLN8 is a positive regulator for the avian adipogenesis.
CLN8 is involved in the gene regulation network of preadipocyte differentiation
To better understand the effect of CLN8 on adipogenesis, we performed mRNA-Seq experiments in CLN8OE, and CLN8NC cells prior to differentiation (day 0) and day 3 after differentiation (Additional file 8, 9). To confirm the results from mRNA-Seq, RT-PCR was performed in the present study firstly. Our RT-PCR results were highly consistent with the mRNA-Seq results (Additional file 10–12). In total, 2488 and 2159 DEGs were up-regulated and down-regulated in CLN8OE cells compared to CLN8NC cells on day 3, respectively (Additional file 13). GO analysis for DEGs up-regulated in CLN8OE cells are enriched for cell differentiation processes, while DEGs up-regulated in CLN8NC cells are enriched for the terms related to cell cycle (Table 2; Additional file 14). These results indicate that the knock-in of CLN8 alone has initiated the preadipocyte differentiation early. Notably, we focused on the DEGs enriched in entries directly related to fat cell differentiation in CLN8OE cells compared to the CLN8NC group. There are 9 up-regulated DEGs were enriched in the positive regulation of fat cell differentiation (PTGS2, PPARG, MAPK14, ASXL2, PPARD, STK4, KLF5, BMP2, ZC3H12A), and 9 down-regulated DEGs were enriched in the negative regulation of fat cell differentiation (ARNTL, RUNX1T1, GATA2, CCN4, YAP1, ANKRD26, ZFPM1, BBS12, ASXL1) (Additional file 14). These data suggest that CLN8 is involved in the gene regulation network of avian preadipocyte differentiation and positively regulates avian adipogenesis by modulating the expression of a set of genes related to fat cell differentiation, including PPARG, which is a core regulator of adipocyte differentiation.
Table 2
The representative GO enrichment analysis terms of DEGs in CLN8OE and CLN8NC cells
Group
|
Up-regulated
|
Gene Count
|
Log(P -value)
|
Down-regulated
|
Gene Count
|
Log(p-value)
|
CLN8OE VS CLN8NC
(Day 0)
|
regulation of cell growth
|
49
|
-11
|
mRNA metabolic process
|
98
|
-17
|
regulation of cell development
|
52
|
-9.8
|
regulation of cell cycle process
|
105
|
-15
|
positive regulation of cellular component biogenesis
|
46
|
-7.4
|
mitotic cell cycle
|
89
|
-14
|
lipid biosynthetic process
|
43
|
-4.6
|
DNA metabolic process
|
102
|
-13
|
positive regulation of cell morphogenesis involved in differentiation
|
11
|
-2.7
|
positive regulation of cell cycle
|
53
|
-7.8
|
CLN8OE VS CLN8NC
(Day 3)
|
positive regulation of cell migration
|
94
|
-19
|
skeletal system development
|
72
|
-13
|
regulation of cell adhesion
|
107
|
-16
|
DNA metabolic process
|
83
|
-9.4
|
positive regulation of MAPK cascade
|
70
|
-11
|
nuclear division
|
46
|
-8.7
|
white fat cell differentiation
|
5
|
-1.9
|
mitotic cell cycle
|
68
|
-8.2
|
positive regulation of fat cell differentiation
|
9
|
-1.2
|
negative regulation of fat cell differentiation
|
9
|
-1.5
|