WGBS of embryonic gonads revealed that long non-coding RNAs in the MHM region might be involved in cell autonomous sex identity and female gonadal development in chickens

Figures & data

Figure 1. DNA methylome characteristics of chicken embryonic gonads: (a) CG methylation level Pearson correlation analysis. The varying shades of color represent the degree of correlation; (b) Principal component analysis (PCA) of DNA methylation at the genome-wide level; (c) genome-wide differences in methylation levels between sexes at E6 (left panel) and E10 (right panel). From outside to inside are male methylation level, difference in methylation level between males and females, and female methylation level, respectively; color grading indicates the methylation level; (d) dendrogram clustering of CG methylation levels in each sample. The closer distance, the higher similarity among them.

Figure 2. Distribution of DNA methylation levels in chicken embryonic gonads: (a) the violin plot displays the DNA methylation levels of CG, CHG, and CHH contexts in males and females at E6 and E10, where H represents any nucleotide other than G; (b) the bar chart illustrates the methylation levels of functional elements, with the gray bars representing T-test error bars. There were no significant differences observed among the samples within each group. CGI (CpG island) refers to DNA regions that have a high density of CG sites and are typically longer (> 200 bp) than surrounding regions. CGI shore represents the regions adjacent to CGIs, generally defined as a DNA region within 1000 bp of the CGI itself.

Figure 3. Cluster analysis of differentially methylated regions between sexes at E6 and E10. The varying shades of color represent the methylation level.

Figure 4. Functional enrichment analysis of differentially methylated genes: (a) the cluster heatmap displays the methylation levels of CG-DMRs between sexes at E6 (left panel) and E10 (right panel). The varying shades of color represent the methylation level; (b) GO enrichment analysis of differentially methylated genes between sexes at E6 (left panel) and E10 (right panel); (c) KEGG pathway analysis of differentially methylated genes between sexes at E6 (left panel) and E10 (right panel).

Figure 5. Functional enrichment analysis of differentially expressed genes: (a) the volcano plot displays DEGs between sexes at E6 (left panel) and E10 (right panel), highlighting the major genes involved in sex determination in chickens; (b) GO enrichment analysis of DEGs between sexes at E6 (left panel) and E10 (right panel); (c) KEGG pathway analysis of DEGs between sexes at E6 (left panel) and E10 (right panel).

Figure 6. Combined analysis of differentially methylated genes and differentially expressed genes: (a-c) the scatter plot shows the log2 fold change of gene expression and differential methylation levels between sexes at E6 (left panel) and E10 (right panel) in promoter (a), exon (b), and intron (c); (d) correlation heatmap of between sexes at E6 (top panel) and E10 (bottom panel) in promoter, exon, and intron; (e) the Venn diagram illustrates the intersection and differences between DMGs and DEGs.

Table 1. Correlation of DNA methylation differences in functional regions with differential mRNA expression.

Stage	Region	Total number of genes	Number of negatively correlated genes	Proportion of negatively correlated genes to the total number of genes.
E6	promoter	32	13	37.5%
E6	intron	357	218	61.06%
E6	exon	38	23	60.52%
E10	promoter	170	85	50%
E10	intron	862	443	51.39%
E10	exon	161	86	53.41%

Table 2. Differentially methylated promoter and differentially expressed genes.

Stage	Region ID	Chr	Difference of methylation (Male vs Female)	Log2Fold change (Male vs Female)	P.value
E6	ENSGALG00000047196	4	11%	4.62	5.54E^−08
	ENSGALG00000001109	27	13%	−1.04	1.59E^−14
	ENSGALG00000051419	Z	54%	−10.92	2.74E^−14
	ENSGALG00000050012	Z	64%	−7.01	3.02E^−05
E10	ENSGALG00000037939	11	25%	−1.85	2.12E^−31
	ENSGALG00000054856	4	20%	−1.87	1.58E^−22
	ENSGALG00000050400	33	21%	−1.08	6.28E^−14
	ENSGALG00000027090	4	20%	−6.54	2.55E^−37
	ENSGALG00000051419	Z	70%	−9.73	1.25E^−15
	ENSGALG00000050012	Z	72%	−5.78	0.00018
	ENSGALG00000048026	1	29%	−3.33	3.68E^−06
	ENSGALG00000052052	6	28%	−1.49	4.22E^−11
	ENSGALG00000050315	1	10%	1.84	2.12E^−05
	ENSGALG00000009452	5	20%	−1.86	3.87E^−11

Figure 7. BSP and qRT-PCR of candidate lncRnas. (a), (b), and (c) show the validation results for ENSGALG00000050012, ENSGALG00000051419, and ENSGALG00000045887, respectively. In each figure, the upper left figure shows the WGBS sequencing result, the box is the verification region, the upper right figure is the BSP sequencing result, the black dots indicate the sites where methylation groups occurred, and the white dots indicate non-methylated sites. The lower figure is the result of qRT-PCR. G, B, and M represent the gonads, brain, and muscle, respectively, and AI represents the aromatase inhibitor treatment group. Asterisks above the bars indicate statistically significant differences between samples: *P < 0.05, **P < 0.01, ***P < 0.005, ****P < 0.001.

Figure 8. FISH reveals the spatial expression of candidate lncRnas. Marker genes and candidate genes are displayed in different colors, with the marker gene CVH for primordial germ cells shown in green, the marker genes FOXL2 for ovarian follicle cells and SOX9 for testicular Sertoli cells depicted in yellow. The candidate genes ENSGALG00000050012 and ENSGALG00000051419 are respectively represented in red and blue. (a) and (b) represent the female gonad and male gonad, respectively. (c) and (d) show enlarged views of the orange-boxed regions in the upper images.

Figure 9. Validation of candidate gene expression after 5-Aza-CdR treatment in DF-1 cells by qPCR. (a): ENSGALG00000050012; (b): ENSGALG00000051419; (c): ENSGALG00000045887. Different letters of a, b and c above the bars indicate a statistically significant difference (p < 0.05) among different groups, while the same letters indicate that the difference is not statistically significant (p > 0.05).

Figure 10. Fluorescence activity analysis of truncated promoter and prediction of transcription factor binding sites. (a) ENSGALG00000050012, (b) ENSGALG00000051419, the top panel and bottom panel were fluorescence activity in promoter and predicted transcription factor binding sites in the promoter, respectively.