m⁶A mRNA Methylation Analysis Provides Novel Insights into Pigmentation in Sheep Skin

Figures & data

Table 1. Primer sequences for RT-Qpcr and MeRIP-Qpcr..

Genes	Primer sequence (5’−3’)	Product length (bp)	Accession NO.
CSPG4	F: CTGGTCCGGCACAAGAAGAT R: AGAACACAATGTCCGCTGGT	109	XM_027957164.2
MAML2	F: GGTTTAACCTCGCCACTCCA R: CCCTTACTTCGGACACTGGT	105	XM_015100757.3
DENND2B	F: CCAGAGCCTCATGGTTCCAG R: TCTTGTTGGCAGTCATGGTCA	108	XM_042233083.1
ADAMTS1	F: CGACAAATGTGGCGTCTGTG R: TTTGTGGCTCCGGTTGGAAT	117	XM_004002797.5
MAP1B	F: CTACGTGGTGAGTGGGAACG R: ATGAGTCGGGATCAGGGTCA	121	XM_027980019.2
GAB2	F: AGAGACAGTGCCTACGACCT R: GCTGGGCGTCTTGAAGGTAT	114	XM_042237732.1
TYRP1	F: TGGCCAGGTGAGTACTGAAA R: CAGAATGGGGTCCCGACAAA	190	EU760771.1
SLC45A2	F: CTCTGGCCATGTGCACCTTA R: GAGAGCCACAAAGCAACAGC	296	XM_004017064.4
TYR	F: GCACAACCGGGAATCCTACA R: CCAGCACAGCAGTAAGGACA	221	NM_001130027.1
TCL1A	F: CCAACCCTGTGTGATCTGCT R: GTATGAGGACCCCGAAGCTG	120	XM_012099174.4
ZAR1	F:GGCACTAACAAGAATTGTAAACAGA R: TTGCCAAACAGCCTTTCACG	303	XM_042251401.1
MCP-3	F: CACACCATCACGGACCAAGA R: CACTGCACATCATCCCTGGT	202	NM_001009411.2
β-actin	F: GCAGGAGTACGATGAGTCCG R: AACCGACTGCTGTCCCCTT	238	DQ152927.1

Table 2. Sequence statistics and quality control.

Sample Name	Raw Reads	Clean Reads	Clean Reads (%)	Q20%)	Q30%)	GC Content (%)
BLACK1-IP	56239292	53025640	88.98	97.47	92.58	53.8
BLACK1-input	65348920	55488324	75.71	97.82	93.51	57.61
BLACK2-IP	57702866	53712660	87.41	97.15	91.75	52.29
BLACK2-input	61572416	53491488	76.40	97.87	93.6	54.53
BLACK3-IP	58047254	53967086	87.495	97.11	91.74	51.59
BLACK3-input	66534292	54903570	73.85	97.78	93.37	57.79
WHITE1-IP	56076312	52831282	88.35	97.65	92.99	54.35
WHITE1-input	64703960	55491504	73.12	98.06	94.1	55.79
WHITE2-IP	59653724	56942474	88.93855	97.74	93.28	53.88
WHITE2-input	55037238	50143426	81.96	97.86	93.57	56.9
WHITE3-IP	53302052	50733478	89.88	97.83	93.54	51.93
WHITE3-input	58191410	52316806	81.33	97.97	93.86	57.66

Figure 1. Peak Distribution. (a) m⁶A peaks distributions in black and white skin of sheep among mRNA transcripts. (b) m⁶A peaks distributions in black and white skin of sheep among genes. (c) Distributions of differential m⁶A peaks in per gene. (d) Numbers of m⁶A peaks in the mRNA transcripts of chromosomes in black skin of sheep. (e) Numbers of m⁶A peaks in the mRNA transcripts of chromosomes in white skin of sheep. (f) the top five motif sequences for m⁶A-containing peaks from all samples.

Figure 2. Peak Distribution. (a) Density distributions of m⁶A peaks in different gene functional elements (5'UTR, CDS, and 3'UTR) in each sample. (b) Distributions of m⁶A peaks in different gene functional elements (5'UTR, 3'UTR, 1 st exon, and other exons) in black skin and white skin.

Figure 3. m⁶A peaks in black skin and white skin of sheep. (a) Significantly different m⁶A peaks in black skin vs. white skin. (b) Significant differences in the distributions of m⁶A peaks on sheep chromosomes.

Table 3. Top 10 significantly upregulated and downregulated m⁶A peaks (BLACK vs. WHITE).

Chromosome	ChromStart	ChromEnd	Gene name	P value	Strand	Log2FC	Regulation	Peak region
4	23383792	23384256	ARL4A	1×10^−13	+	−3.94	down	CDS
21	19123523	19125178	TENM4	7.9×10^−21	+	−3.94	down	CDS
18	31064789	31065785	CSPG4	6.3×10^−12	+	−3.70	down	CDS
15	15665020	15666219	MAML2	6.3×10^−16	+	−3.63	down	3‘UTR
15	46957631	46983543	DENND2B	7.9×10^−17	+	−3.62	down	CDS
5	40849569	40849904	TRIM52	1.3×10^−11	-	−3.62	down	CDS
2	2.24E+08	2.24E+08	PLEKHM3	7.9×10^−15	-	−3.43	down	5’ UTR
23	53963974	53964273	CTIF	5×10^−14	+	−3.43	down	3’ UTR
16	4376992	4379027	SH3PXD2B	5×10^−13	-	−3.38	down	3’ UTR
18	68897002	68898246	BAG5	4.9×10^−10	-	−3.36	down	CDS
11	9176017	9189954	TNRC6C	2. ×10^−19	+	3.95	up	CDS
14	52705416	52706353	ITPKC	5×10^−13	-	3.86	up	-
21	42866241	42867012	UQCC3	6.3×10^−14	+	3.72	up	3’ UTR
13	25201358	25202258	OTUD1	2×10^−21	+	3.67	up	CDS
11	27286934	27287133	ANKRD40	2.5×10^−11	-	3.62	up	3’ UTR
14	54169492	54170493	CIC	1.4×10^−10	-	3.60	up	CDS
21	10853343	10854221	CREBZF	1.1×10^−8	+	3.57	up	CDS
18	47442175	47443157	CLEC14A	1.7×10^−9	-	3.53	up	CDS
16	74060446	74061095	ICE1	1×10^−15	-	3.51	up	CDS
3	5933126	5933674	ABL1	3.2×10^−12	-	3.49	up	CDS

Figure 4. GO analyses of m⁶A modified genes. Top 15 GO terms of the differentially methylated downregulated genes and upregulated genes in three categories (BP, MF, CC). The red column indicates significant enrichment, while the blue column indicates insignificant enrichment.

Figure 5. KEGG analyses of m⁶A modified genes. The enriched KEGG pathways of the differentially methylated downregulated genes (a) and upregulated genes (b).

Figure 6. GO and KEGG analyses of DEGs. (a) Volcano plot showing the differential gene expression in black skin vs. white skin of sheep. (b) Top 20 KEGG pathways enriched for DEGs. (c) Top 15 GO terms of DEGs in three categories (BP, MF, CC). The red column indicates significant enrichment, while the blue column indicates insignificant enrichment.

Table 4. Top 10 significantly upregulated and downregulated genes (BLACK vs. WHITE).

Locus	Gene name	P value	Strand	Log2FC	Regulation
chr18:60855163–60858729	TCL1A	0.034	-	−4.55759	down
chr 6:73764050–73767864	ZAR1	0.026	+	−4.5438	down
chr 18:32645618–32716451	MCP-3	6.37×10^−7	+	−4.12923	down
chr 7:75675023–75677892	SIX1	0.008	-	−4.117	down
-	5S_rRNA	0.022	-	−2.73541	down
chr 1:97071826–97207889	SYCP1	0.019	+	−2.73325	down
chr 1:219010098–219035477	KNG1	0.019	-	−2.53527	down
chr 3:146247496–146250747	AQP5	0.016	-	−2.45983	down
chr 14:56880795–56886640	FOXA3	0.036	+	−2.2333	down
chr 3:39631341–39658950	CAPN14	0.033	-	−2.20124	down
chr 2:87540480–87556594	TYRP1	6.9×10^−14	+	13.08397	up
chr 16:41978919–42013172	SLC45A2	3.57×10^−11	+	9.966668	up
chr 11:11951667–11954992	TSPAN10	1.68×10^−10	+	9.599748	up
chr 2:79584125–79594791	MLANA	7.04×10^−12	+	8.412334	up
chr 18:25349881–25450837	TRPM1	1.91×10^−7	-	7.728241	up
chr 3:175255140–175263368	PMEL	2.16×10^−16	+	6.455803	up
chr 21:7263507–7379923	TYR	2.33×10^−11	-	5.898827	up
chr 5:5460092–5524649	UNC13A	0.004	+	5.244571	up
chr 8:12050576–12225332	SNAP91	0.0008	-	5.193136	up
chr 1:8353926–8356212	KCNJ13	0.006	+	4.950321	up

Figure 7. Combined m⁶A-seq and RNA-seq analysis. (a) the overlap gene between downregulated m⁶A methylation and downregulated mRNA in black skin vs. white skin (b) Four-quadrant diagram depicting the distributions of m⁶A modified genes and DEGs. (c-d) the enriched KEGG pathways among the identified (c) hyper-up genes and (d) hypo-up genes.

Figure 8. Verification of DEGs and DMGs. (a) the m⁶A methylation modifications of six genes verified by MeRIP-Qpcr.. (b) the expression levels of six genes verified by qPCR.

Data availability statement

All datasets used in this study are available from the corresponding author on reasonable request.