Differences in the fast muscle methylome provide insight into sex-specific epigenetic regulation of growth in Nile tilapia during early stages of domestication

ABSTRACT

Growth is a complex trait whose variability within a population cannot be explained solely by genetic variation. Epigenetic regulation is often suggested as an important factor shaping the phenotype, but its association with growth can be highly context- and species-dependent. Nevertheless, the mechanisms involved in epigenetic regulation of growth in fish are poorly understood. We have used reduced representation bisulphite sequencing to determine the genome-wide CpG methylation patterns in male and female Nile tilapia of different sizes but at the same early stage of domestication. The average CpG methylation level in the reduced genome representation was 63% across groups but many sites displayed group-specific methylation patterns. The number of differentially methylated (DM) CpGs was much higher when the interaction between sex and weight was included rather than when these factors were considered separately. There were 1128 DM CpGs between large and small females and 970 DM CpGs between large and small males. We have found many growth-related genes associated with DM CpGs, namely map3k5 and akt3 in females and gadd45g and ppargc1a in males. Only 5% of CpG locations associated with growth were common to both sexes. In particular, the autophagy-related gene atg14 displayed a high association of methylation with growth exclusively in males. The sexually dimorphic association between atg14 methylation and growth may uncover novel metabolic mechanisms at play during mouth brooding in Nile tilapia females. Taken together, our data suggest that epigenetic regulation of growth in Nile tilapia involves different gene networks in males and females.

KEYWORDS:

• DNA methylation
• RRBS
• sexual dimorphism
• teleosts
• muscle growth
• Oreochromis niloticus

Background

Growth in fish is a polygenic and complex trait, which involves a panoply of pathways and tissues [1] and is influenced by the interaction with other fitness-related traits [2]. The role of abiotic and biotic factors such as temperature, water current, density, diet, and social interactions in shaping fish growth or condition factor is well established, with ample experimental evidence from wild and farmed fish populations [3–5]. In wild populations, in particular, stochastic conditions such as food availability or disease outbreaks and additional selective pressure (e.g., predation) can lead to expression of sub-optimal growth rates or prevent accurate predictions of growth [6]. In contrast, captive populations benefit from relaxed natural selection and therefore can exhibit wider phenotypic variation than their wild conspecifics [7]. First, low mortality in aquaculture conditions allows for survival of phenotypes that would not persist in the wild [8]. Second, specific conditions of the controlled rearing environment (food, density, temperature) are often aimed at promoting growth rate, i.e., allowing individuals with high growth potential a less restricted expression of this trait. Relaxed natural selection, when leading to increased phenotypic variation [9], can help the subsequent artificial selection and/or rapidly occurring adaptation to captivity, which can be seen as one of the first steps in the domestication process [10].

The initial stage of domestication represented by a parental generation (F0) produced in the wild but reared under farmed conditions is a unique model for studying growth mechanisms. It represents a particular case of reaction norms generated by a new and highly controlled environment, where individuals benefit from reduced selective pressure from the wild but are not yet affected by artificial selection. Indeed, man-assisted selection such as breeding of fast growers is often directional and stabilizing, thus domestication-related changes in the subsequent generations will be rather associated with relocation and/or flattening of reaction norms of selected and co-selected traits [11].

Stable changes in phenotypic expression can be extremely fast under captivity conditions, even within a single generation of domestication [12] and do not always require neo-Darwinian modulation of allelic frequencies [13]. Epigenetic changes such as DNA methylation, i.e., potentially heritable modifications of chemical structure of genome without affecting its nucleotide composition, can be an important mechanism shaping phenotypic differences occurring within the same generation. Cytosine methylation is among the most studied epigenetic marks [14] and it has recently been associated with hatchery rearing in coho salmon (Oncorhynchus kisutch) [15]. Global changes in DNA methylation have been observed in Nile tilapia (Oreochromis niloticus) gonads after temperature treatment [16] and it has been shown that differential cytosine methylation may contribute to sex-specific growth phenotypes in hybrid tilapia (O. mossambicus × O. niloticus) [17]. Differences in cytosine methylation of the myogenin gene have also been linked to thermal plasticity of growth in Atlantic salmon (Salmo salar) and Senegalese sole (Solea senegalensis) [18,19]. Therefore, a better knowledge on phenotype-associated methylome modifications may have key implications in fish epigenetic programming in the future [20].

In this study, we assessed the size-related differences in genome-wide methylation status between F0 generation male and female Nile tilapia (O. niloticus), which were born in the wild and reared under hatchery conditions. Growth of this important aquaculture species is a highly plastic [21] and heritable trait [22]. Thus, detecting methylation changes related to extreme growth levels can undercover new regulatory mechanisms of myogenesis and provide potential epigenetic biomarkers of such trait. Moreover, due to the magnitude of sexual size dimorphism in Nile tilapia (with males being considerably larger than females) [23], we hypothesized that there would be sex-specific differences between methylation profiles in both large and small individuals. By using the wild-born but hatchery-reared individuals, we are simulating the first steps of domestication, i.e., acclimation to farmed conditions and expression of reaction norms that are not necessarily attained in wild conditions. Such a less restricted expression of growth potential allows us to get insight into potential epigenetically regulated functions specific to growth and to sex in fish.

Results

Experimental groups and library characterization

Within full-sib groups reared in common garden, the large individuals were 8.7- and 5.7-fold heavier than their small counterparts, in the case of males and females, respectively (Figure 1). A mean of 30 M raw reads was obtained per library, of which 15 M were uniquely mapped (i.e., 50% mapping efficiency) and an additional 8 M reads could be multiple mapped (Figure 2, Supplementary Table 1). Neither individual samples nor compared groups differ significantly with respect to their raw, trimmed, or mapped reads. The mean CpG coverage of reads mapping uniquely to the fraction of the genome that was successfully aligned with reads is 31-fold or higher (Supplementary Figure 2). Absence of a peak at the right side of the histogram indicates that the data did not contain polymerase chain reaction (PCR) duplicates.

Figure 1. Average weight of Nile tilapia, according to sex and weight (n = 6). Large females were statistically heavier than small females (t-test, p-value <0.001) and lighter than large males (t-test, p-value <0.001), while small males were statistically lighter than large males (t-test, p-value <0.001) and heavier than small females (t-test, p-value <0.001).

Figure 2. Number of total raw (dark grey), quality-trimmed (light grey), adapter-trimmed (purple), uniquely mapping (green), and multiple mapping (red) reads in the group of large and small females and males, respectively.

The fast muscle methylome in Nile tilapia

The global cytosine count corresponded to 244 M reads with a fold coverage was 3.2. Average methylation corresponded to 63% (Table 1), mean methylation levels per chromosome were homogeneous (Supplementary Table 2), and the distribution of methylated cytosines was similar across all samples (Figure 3). However, methylation levels differed between chromosomes, namely between chromosome LG3b and LG4 (Figure 3). In addition, the genomic context of the RRBS dataset differed from that in the whole Nile tilapia genome (Figure 4). Indeed, promoter-annotated cytosines were enriched twofold in the reduced representation genome compared to whole genome. A twofold increase in CpG sites found in the RRBS dataset compared to the whole genome was also observed for exon regions.

Table 1. Global cytosine methylation level in CpG, CHG, and CHH contexts; number of cytosines analysed and cytosine coverage with respect to the whole Nile tilapia genome.

Weight	Sex	Cytosine count (10^6)	Cytosine fold coverage	% mCpG	% mCHG	% mCHH
Large	Female	243.3 (±12.7)	3.2	64.1	1.3	1.1
	Male	254.2 (±19.9)	3.3	60.6	1.2	1.0
Small	Female	241.4 (±21.4)	3.1	64.4	1.3	1.0
	Male	236.1 (±15.4)	3.0	62.1	1.3	1.0

Figure 3. Heat map of methylation levels of 234,755 CpG sites distributed across the Nile tilapia linkage groups (columns) and separated in sex- and weight-specific groups (rows, n = 6). Methylation levels below 40%, between 40% and 60%, and above 60% are indicated in blue, yellow, and red, respectively.

Figure 4. Genomic context of CpG sites found in the whole genome (a) and in the reduced representation dataset (b).

Sex- and weight-specific differentially methylated (DM) cytosines

Unequal distribution of DM CpG sites between sexes

A total of 234,755 CpG positions that were covered in every sample at least once were used for analysis of differential methylation between sexes. Interestingly, neither genomic context (Figure 5) nor chromosomal distribution (Supplementary Table 3, Supplementary Figure 3) were similar between the sets of hyper- and hypomethylated CpG sites. For example, transcription termination sites were represented by 10% of hypomethylated CpG sites in males but absent in hypermethylated CpG sites.

Figure 5. Genomic context of sex-specific differentially methylated cytosines (hypo- and hypermethylated). Twelve fish per sex and 234,755 cytosine positions in total were compared.

Moreover, compared to the CpG sites available in the RRBS dataset (Figure 4) there was an approximately twofold increase in the proportions of all sex-related DM CpG sites (Figure 5) found in promoters and a sixfold decrease in those found in exons. In general, the distribution of DM CpG sites across the genome was very heterogeneous (Supplementary Table 3). Table 2 shows the genes that are associated with DM CpGs and that are functionally relevant to muscle-related or epigenetics-related processes. Out of a set of selected genes connected to DM CpG between males and females, the gene encoding for atg14 (autophagy-related protein 14) is linked with four exonic CpG sites, all hypomethylated in males. Other hypomethylated DM CpG sites in males were associated with hdac11 (histone deacetylase 11), tgfb2 (transforming growth factor beta 2), and lepr (leptin receptor) and one hypermethylated CpG site is connected with itga4 (integrin subunit alpha 4).

Table 2. List of functionally relevant genes associated with CpG sites that are differentially methylated between sexes in Nile tilapia.

Chr	Start	Diff. (%)	Strand	Annotation	Distance nearest TSS	RefSeq Nearest promoter	Gene symbol	Function	DM CpG count	Frequency ratio
LG17	36204426	−25.10	+	TTS	17,789	XM_019346989.1	lepr	Myoblast differentiation, fat deposition	1	0.13
LG16	15827762	27.94	+	Intergenic	−15,883	XM_019353317.1	itga4	Myogenesis	1	0.11
LG5	2724018	−36.77	-	Intergenic	−2827	XM_003446517.4	hdac11	Histone modification, myoblast differentiation	1	0.14
LG19	1368191	−31.59	+	Exon	−2614	XM_003457198.4	atg14 autophagy-related		4	0.17
LG19	1368164	−29.91	+	Exon	−2641				4	0.17
LG19	1368175	−29	+	Exon	−2630				4	0.17
LG19	1368168	−26.83	+	Exon	−2637				4	0.17
LG11	29219134	−33.50	+	Intron	65,818	NM_001311314.1	tgfb2	Muscle skeletal atrophy	2	0.06

Twelve males were compared to 12 females. Only the results with q-value <0.01 are retained. Positive percentage values indicate hypermethylation in males when compared to females, whereas negative values indicate hypomethylation in the same comparison. Frequency ratio corresponds to a number of unique gene-associated DM CpG sites divided by a total number of CpG sites available for the analysis and associated with the same gene.

Substantial differences in CpG methylation between large and small individuals of the same sex

Out of approximately 300k CpG locations, a total of 1128 and 970 DM CpG positions (false discovery rate (FDR) <1%, difference >25%) were found in large individuals when compared to small ones within females and males, respectively (Supplementary Table 4a). The overall genomic context in both sets of DM CpG positions was similar (Supplementary Figure 4). However, the distribution of DM CpG sites across chromosomes differed between male- and female-specific comparisons (Figure 6). Chromosomes 2, 15, and 16 are the most enriched in DM CpGs in females, and chromosomes 2, 3a, 15, and 22 in males (Figure 6). In males, chromosomes with the highest number of DM CpG sites also contained the sites with the highest values of differential methylation, particularly hypermethylated sites. The shortest linkage group (LG3a) contains one of the highest proportions of DM CpG per nucleotide, whereas LG3b shows one of the lowest proportions of DM CpG per nucleotide.

Figure 6. Chromosomal distribution of weight-specific differentially methylated CpG sites in females (1128 positions) and males (970 positions). Small individuals are set as reference for both sexes separately. Methylation differences in males are represented in the inner circle (blue background) and those of females in outer circle (grey background). Histograms pointing inwards (dark green) and outwards (dark red) represent hypomethylated and hypermethylated sites, respectively.

A list of genes associated with growth-related DM CpG loci in females (Supplementary Table 4a) and males (Supplementary Table 4b) was created based on their functional relevance to muscle- or epigenetics-related processes. Among the selected genes associated with DM CpG in females, mettl9 (methyltransferase-like protein precursor 9) is associated with the highest number of DM CpG (five in total for variant X4 and X8 of mettl9), all hypomethylated in large females and close to the transcription start site (TSS) of that gene. Their frequency ratio (number of DM CpG linked with the specific gene/number of all CpG sites related to the gene and available in the RRBS data) was higher than 50%. A gene encoding for ppp1r3c (protein phosphatase 1 regulatory subunit 3C) was represented by 2 DM CpG sites, one hypomethylated and situated in exon and another one hypermethylated and situated in intergenic region. Similarly, histone deacetylase gene (hdac11) was also represented by two DM CpG, one hypo- and one hypermethylated. We observed hypermethylation in exonic regions of gata6 and xpo5 (exportin 5), in intergenic regions close to dmap1 (DNA methyltransferase 1-associated protein) and intronic region of ranbp2, encoding for a SUMO E3 ligase (RAN binding protein 2). Two genes, map3k5 (mitogen-activated kinase kinase kinase 5, known also as ask1) and esrrg (encoding for oestrogen-related receptor gamma), were affected by hypermethylation of multiple CpG sites (five and four, respectively) in their promoter regions, very close to their TSS. In contrast, CpG hypomethylation was detected in introns of igf2bp2 (insulin-like growth factor 2 mRNA binding protein 2 (IGF2BP2)), dner (delta/notch like EGF repeat containing gene), and itgbl2 (integrin beta-like 2), as well as in the intergenic region close to prlr (prolactin receptor).

In males, the highest number of DM CpG was associated with atg14 (autophagy-related gene 14). This gene was represented by 20 DM CpG (87% of all CpG sites associated with atg14 and available in the RRBS dataset), all hypomethylated in large males and situated in introns or exons of tbp2 (TATA-box binding protein 2). Other hypomethylated CpG sites were associated with ppargc1a (peroxisome proliferator-activated receptor-γ coactivator-1α, also called PGC-1alpha), in the promoter region of gdf3 (growth differentiation factor 3) and in the third exon of gadd45a (growth arrest and DNA damage-inducible protein GADD45 gamma-like). We also observed hypomethylation of five CpG sites situated in the promoter region of hadh (hydroxyacyl-CoA dehydrogenase) and two CpGs located in the exon region of phf20 gene (PHD finger protein 20), as well as hypermethylation of CpG sites associated with gata6 and xpo5.

Only a few growth-related DM CpG sites are shared between males and females

Out of approximately 1000 possible DM CpG locations, only 57 (5.7%) were common to both sexes. Among them, 28 showed the same methylation changes in both sexes (i.e., if hypermethylated in big females, they were also hypermethylated in big males) (Supplementary Figure 5, Supplementary Table 4b). Surprisingly, the remaining (29) common DM CpG locations showed opposite patterns of methylation associated with growth (Supplementary Table 4b). Several genes common for female- and male-specific analysis were selected based on their functional relevance to muscle- or epigenetics-related processes (Table 3) and showing similar or opposite patterns of methylation related to growth. Among them, the aforementioned gata6 was related to hypermethylated DM CpGs in both large males and large females. In contrast, xpo5, a gene encoding for Exportin 5, was associated with a CpG site significantly hypermethylated in large males, but hypomethylated in large females.

Table 3. List of functionally relevant genes associated with weight-related DM CpG sites that are found separately in both sexes in Nile tilapia.

Chr	Start	Diff. females (%)	Diff. males (%)	Strand	Annotation	Distance nearest TSS	RefSeq Nearest promoter	Gene symbol	Function
LG18	25283313	25.60	28.36	−	Exon	6115	XM_005476460.3	gata6	Transcription factor, smooth muscle growth
LG17	29592811	25.79	28.02	+	Intron	−1607	XM_019346583.1	serpinb	Enriched in dystrophic muscle
LG6	3091542	25.34	−27.15	−	Intron	−17,945	XM_019359489.1	sorcs1	Energy balance in brain
LG14	16673925	−29.84	−27.86	−	Exon	1239	NM_001311336.1	neu4	Sialic acid removal
LG13	9351833	−34.17	26.59	+	Exon	5364	XM_003451992.3	xpo5	Export of miRNA form nucleus to cytoplasm

Only the results with q-value <0.01 are retained. Positive percentage values indicate hypermethylation, whereas negative values indicate hypomethylation. For each sex, the small group was used as reference.

The comparison between large and small fish yields few DM CpGs when sex information is omitted

A total of 117 DM CpG sites were found when the 12 largest individuals were compared to the 12 smallest individuals from both sexes (Supplementary Table 4a). Their distribution was uneven and different from sex-specific DM CpG distribution (Supplementary Figure 6).

Table 4 shows a subset of weight-specific CpG sites associated with genes involved in growth in fish. Gata6 contained the highest number of DM CpG sites (4), all exonic and hypermethylated. A gene encoding for igf2bp2 (IGF2BP2) was associated with one intronic hypermethylated DM CpG, and itgbl2 was associated with one intergenic hypomethylated CpG. The first gene was also identified in the size-related comparison within females, whereas the second was found both within females (Supplementary Table 5a) and males (Supplementary Table 5b). Only a few unique genes were associated with DM CpG sites in this group, namely cdh2 (cadherin 2) and hadh, a gene that was previously found to be size-related in males.

Table 4. List of functionally relevant genes associated with DM CpG sites related to growth in Nile tilapia.

Chr	Start	Diff. (%)	Strand	Annotation	Distance nearest TSS	RefSeq Nearest promoter	Gene symbol	Function	DM CpG count	Frequency ratio
LG17	29592812	26.78	+	Intron	−1607	XM_019346583.1	serpinb1	Enriched in dystrophic muscle	1	0.17
LG18	25283253	27.55	−	Exon	6054	XM_005476460.3	gata6	Transcription factor, smooth muscle growth	4	0.12
LG18	25283314	26.95	−	Exon	6115				4	0.12
LG18	25283294	26.64	−	Exon	6095				4	0.12
LG18	25283290	26.67	−	Exon	6091				4	0.12
LG16	25397159	34.60	+	Intron	−2578	XM_019352870.1	igf2bp2	Muscle differentiation	1	0.10
LG15	2856820	−28.34	+	Intergenic	−5997	NM_001279735.1	thbs1a	Insulin resistance	1	0.05
LG13	30916611	33.35	+	Intergenic	−12,726	XM_019366349.1	hnrnpll	Regulates stat5a involved in growth	1	0.04
LG16	2999356	−25.63	+	Intergenic	−13,226	XM_019353331.1	itgbl2	Hypertrophy, catabolism	1	0.02
LG11	10465297	25.46	+	Exon	68,928	XM_013264857.2	cdh2	Muscle differentiation	1	0.01
LG6	22757890	−26.07	+	Promoter-TSS	−18	XM_003452225.4	hadh	Autophagy, muscle degradation	1	0.09

Twelve large individuals were compared to 12 small individuals. Only the results with q-value <0.01 are retained. Positive percentage values indicate hypermethylation in large individuals when compared to small individuals, whereas negative values indicate hypomethylation in the same comparison. Frequency ratio corresponds to a number of unique gene-associated DM CpG sites divided by a total number of CpG sites available for the analysis and associated with the same gene.

A large number of DM CpG-associated genes are related to growth or involved in epigenetic regulation

Although Functional Annotation analysis did not yield any significant enrichment (Supplementary Table 6), many manually annotated genes were functionally related to growth and epigenetic processes (Figure 7). In the comparison between males and females, we highlighted hdac11, a gene coding for a key epigenetic histone modification factor (histone deacetylase) that is capable of decreasing myoD-dependent transcription [24]. It is important to notice that histone deacetylases were also found in other comparisons. Among other sex-related and DM CpG-associated genes, lepr can have a role in proliferation and differentiation of skeletal myoblasts [25] whereas itga4 is involved in myogenesis [26]. Muscle atrophy was a common function of tgfb2, a transforming growth factor [27], and of atg14, from a family of autophagy-related genes involved in muscle degradation [28,29]. Interestingly, the last gene was also associated with many DM CpGs between large and small males, but not between large and small females.

Figure 7. Schematic representation of genes associated with differentially methylated CpGs in the comparison between males (a), females (b), growth (c), or sex (d) and that play a role in muscle growth or epigenetic regulation processes. Genes in black are specific to one comparison, whereas genes in red are associated with differentially methylated CpGs in more than one comparison. In the male-specific comparison, gadd45g is connected to several other differentially methylated genes. Gadd45g activates the p38MAPK pathway, which is important in muscle growth and regeneration.

In females, a comparison between large and small fish showed DM CpG-associated genes whose function was related to growth, myogenesis, or mitochondrial activity (Supplementary Table 4a). For example, ppp1r3c is involved in glycogen metabolism [30], itgbl2 (Integrin subunit beta-like 2) might contribute to muscle hypertrophy [31], dner (Delta and Notch-like epidermal growth factor-related receptor) plays a role in signalling in muscle stem cells [32], neo1 (encoding Neogenin) promotes myotube formation and regulates myofiber size [33], map3k5 (Mitogen Activated Protein Kinase Kinase Kinase 5) promotes myogenic differentiation [34], esrrg regulates mitochondrial activity [35], and igf2bp2 (Insulin-like growth factor 2 binding protein 2) is involved in myoblast proliferation and myogenesis [36]. The latter can also be involved in epigenetic regulation [37]. Epigenetic processes were also highlighted by other genes, such as xpo5 involved in transfer of miRNA from nucleus to cytoplasm [38] and dmap1 involved in global maintenance methylation [39]. In addition, ranbp2 can modify histone deacetylase HDAC4 [40].

A male-specific comparison between large and small fish showed similar genes related to epigenetic regulation (Supplementary Table 4b), such as hdac9, xpo5, and phf20, a gene encoding for PHD finger protein 20 and involved in histone acetylation [41]. Many DM CpG sites were associated with atg14, but situated in the last introns and exons of the vertebrate-specific tbp2 gene [42] (Figure 8) encoding for TATA-box binding protein 2, required for transcription initiation of promoters that contain a TATA-box [42]. Other growth-related and DM CpG-associated genes found in the male-specific comparison were hadh involved in catabolism, gadd45g (growth arrest and DNA damage-inducible protein gamma-like) causing muscle atrophy [43], and gdf3 and ppargc1a (encoding for peroxisome proliferator-activated receptor gamma) both involved in muscle metabolism [44,45].

Figure 8. Location at 10k (a) and 1k (b) base pair resolution of 20 DM CpG sites situated within introns and exons of tbp2 (TATA-box binding protein 2) and associated with the transcription starts site (TSS) of atg14 (autophagy-related gene 14).

The comparison between large and small fish regardless of sex (Table 4) revealed a growth-related gene cdh2 (encoding cadherin 2) that promotes skeletal muscle differentiation [46] and was only identified in this comparison. In addition, the transcription factor gata6 was found both in size- and sex-specific comparisons between large and small fish (Figure 7), and is involved in regulation of growth in smooth muscle [47].

Discussion

In this study, we report the reduced representation profiling of the fast muscle methylome in Nile tilapia groups with different growth rates. All fish were first generation of captive-reared descendants of wild parents; therefore, their phenotypic expression is on the one hand unrestricted by selective pressures present in wild but absent in hatchery conditions, and on the other hand not yet affected by artificial selection due to selective breeding. To the best of our knowledge, this is the first single-cytosine resolution study of methylation in fast muscle of Nile tilapia in the context of sex, growth, and the early domestication process in fish. Previous studies on genome-wide Nile tilapia methylation were performed on gonadal tissue and at a lower resolution using MeDIP technique, then bisulphite treatment of a restricted number of genes [48]. To date, there is only one report on genome-wide methylation patterns in tilapia [17], suggesting that DNA methylation in skeletal muscle contributes to sex-specific phenotypes. However, this whole genome bisulphite sequencing study focused on a hybrid tilapia species (O. niloticus × O. mossambicus) and the analysis was performed at the level of DM regions, which fails to detect unique but highly DM cytosines and to annotate them separately to specific genes. Importantly, it focused only on sex-specific differences, without taking into account the growth/size of the fish, which is known to be sexually dimorphic in most tilapia species. Therefore, some differences in methylation specifically due to size might have been mistakenly highlighted as sex-specific. In contrast, our study established a multiple factor analysis of methylation data at single nucleotide resolution for sex and growth in combination and independently. Some genes or gene families found to be associated with differential methylation between sexes in hybrid tilapia [17] were also found in Nile tilapia, but not necessarily with the same patterns of methylation. For example, a member of solute carrier family 22 (Supplementary Table 4a) previously found to be hypermethylated in males compared to females was linked in the present study to two hypermethylated CpGs only in large females compared to smaller ones. A gene from the family of histone deacetylases, hdac4, was associated with a hypermethylated region in female hybrid tilapia compared to males [17]. In our study, we detected two histone deacetylases, hdac9 (with hypermethylated CpG sites in large males compared to small males) and hdac11. The latter was associated with hypomethylated CpG in males compared to females, both in hybrid tilapia and in Nile tilapia, but it was also linked to two hypermethylated CpG sites in large females compared to small ones in the present study. It is important to keep in mind that lack of detection of a specific gene by RRBS does not mean that differences in methylation do not exist, since the data represent only a fraction of the whole genome. Nevertheless, our study represents a complementary and more accurate insight into sex- and growth-specific methylation patterns and provides a list of DM CpG-associated genes relevant to muscle growth and epigenetic processes, found specifically or commonly to growth, sex, or interaction between both (Figure 7). As the fish were reared in the same tank, we cannot exclude the possibility that differences in growth rate may be attributed to nutrient condition and feeding behaviour, as well as epigenetic changes.

Our results indicate that global methylation levels in Nile tilapia fast muscle are above 60%, similar to those found in brain of zebrafish [49] of skeletal muscle of tilapia [17]. Moreover, the alignment efficiency of raw reads was on average 50% (Figure 2, Supplementary Table 1), which is in accordance with another study using the same protocol and alignment method in zebrafish [49]. The optimized double-digestion step of our protocol allowed us to study a larger portion of the genome and the choice of Taqα1 in addition to MspI extended the coverage of high CpG density regions as shown in mammalian genomes [50] and in our twofold enrichment of promoter-annotated cytosines. The methylation levels are highly heterogeneous between the chromosomes, and these non-random methylation patterns can be placed in the genomic context (inter- or intragenomic region, promoter, intron exon, or transcription termination site) and assigned to a specific gene. However, it is difficult to predict how methylated cytosines could affect the gene only based on its genomic context. For example, in contrast to generally repressive role of higher methylation near the TSS, gene body methylation can increase or at least not preclude its expression [51]. The order of exons affected by methylation influences gene expression differently [52] and intragenic DNA methylation can also provide splicing signals [53] and suppress cryptic promoters found within transcription units [54]. Finally, CpG methylation can impact distal genes in 5ʹ and 3ʹ direction or act as distal enhancer of one gene while being promoter of another one [55], meaning that transcriptional information would be needed to predict a regulatory outcome of DNA methylation. Our study lacks the experimentally validated connection between changes in methylation of specific DM CpG sites and their consequence in terms of gene expression or phenotype. However, the known functions of DM genes give us an insight into the potential consequences on growth or growth-related processes. As automated functional annotation analysis did not yield any significantly enriched biological processes associated with genes containing DM CpG sites (Supplementary Table 6), we manually selected the genes with growth-related or epigenetics-related functions based on available literature. To our knowledge, none of the available tools for systematic analysis of functional annotation are adapted to RRBS datasets, since they do not consider the number of DM CpG sites associated with a single gene nor the level of methylation differences between compared groups.

Global methylation levels are similar between males and females

The sex-dependent comparison is interesting from an aquaculture point of view, since in Nile tilapia, as in most tilapiine species, males naturally grow faster and attain larger size than females [56]. Factors related to spawning (nutrient loss during spawning, fasting, and post-hatching) have been proposed to explain such dimorphism [57], but epigenetic mechanisms that could affect muscle growth in a sex-dependent manner were never reported. Only 0.01% of all compared CpG sites were DM between males and females and mostly hypomethylated in males. The highest number of DM CpG sites was found in chromosome LG2, which is not the chromosome bearing the highest number of sex-linked markers [56]. Interestingly, while the reduced representation dataset is twofold enriched in promoter- and exon-related CpGs compared to those representing the whole genome, sex-related DM CpGs are even more frequent in promoters but extremely low in exon regions when compared to the raw RRBS dataset. Therefore, methylation upstream of the TSS might be a more frequent sex-dependent mechanism of transcriptional regulation than the methylation within the gene body.

Although functional analysis did not show any significant enrichment of specific functions (Supplementary Table 6), several sex-dependent DM CpGs were associated with genes that have a known function in muscle-related processes (Table 3). Hdac11 (histone deacetylase) is an important epigenetic histone modification factor that can decrease myoD-dependent transcription and thus inhibit myoblast differentiation [24]. Indeed, elevated expression of hdac11 decreased the levels of histone acetylation at promoter regions of mef2c and myoG, which prevented myo-D-dependent myoblast differentiation. Hdac11 is also involved in many immune functions [58]. Lepr can have a role in proliferation and differentiation of skeletal myoblasts [25] and is a potential candidate for growth rate and fat depositions traits in Durocs pigs [59] and is DM between normal and heat-stressed pigs which can affect muscle development and meat quality [60]. Sexual dimorphism of leptin receptors expression in skeletal muscle was also found in human [61]. Itga4 plays a role in myogenesis [26] and is highly downregulated in Myo-deleted mice [62]. Tgfb2 is a transforming growth factor that plays an important role in skeletal muscle atrophy [27]. Four CpG hypomethylated in males were found in exon of atg14, a key component in autophagy, a process of cellular self-renewal [28]. Atg14 also plays an important role in hepatic lipid metabolism [63] and can, by similarity to other autophagy-related genes, have a role in regulating muscle mass due to its potential role in muscle degradation as shown in case of fasting rainbow trout [29].

The importance of sex in the growth-related muscle methylome

Although body weight differences between small and large groups were smaller in females than in males, females displayed more growth-related DM CpG than males. However, their genomic context was very similar between both sexes and close to the genomic context of all CpG sites found in RRBS dataset. In addition, the chromosomal distribution of DM CpGs is strikingly different between sexes, suggesting that different growth-related gene networks are epigenetically modulated in female and male of Nile tilapia. Indeed, out of approximately 1000 possible common positions of sex-related DM CpG sites, only 57 were common to both sexes. For example, the promoters of mettl9 and mettl16 methyltransferase-like genes absent in the DM CpG male dataset were highly represented by hypomethylated CpG in large females. The function of proteins encoded by these methyltransferase-like genes is not well known but mettl16 is suggested to catalyse the methyl marks on RNA, another epigenetic regulator [64] and the mettl3, another gene from methyltransferase-like group can promote muscle cell differentiation by affecting myoD transcription [65]. Large females were also represented by two CpGs associated with ppp1r3c, which is associated with energy metabolism in skeletal muscle [30]. Indeed, ppp1r3c is also called protein targeting to glycogen, since it has an important role in glycogen metabolism [66]. Interestingly, the two DM CpGs showed distinct features. One CpG was hypermethylated and found in an intergenic region whereas the other one, hypomethylated, was found in one exon. In addition, myogenesis, growth, and mitochondrial activity were the main functions of many more genes associated with DM CpG between large and small females, as shown in Figure 7.

Another important gene with two hypermethylated intronic CpG sites in large females encodes IGF2BP2. Its methylation is inversely associated with its mRNA abundance [67]. In addition, IGF2BP2 is involved in satellite-cell-mediated skeletal muscle differentiation, proliferation, and regeneration [36]. Moreover, IGF2BP2 is reported to enhance mRNA stability through recognition of RNA N 6-methyladenosine [37]. Several genes coding for proteins involved in epigenetic regulation were also DM in females (Figure 7). Among them, HDAC11 was discussed above, since it was also associated with sex-related DM CpG sites. Exportin 5 is involved in export of miRNA from nucleus to cytoplasm [38] and DMAP1 is involved in global maintenance methylation [39]. In addition, ranbp2 encodes for a SUMO E3 ligase that is known to modify histone deacetylase HDAC4 [40]. Cytosine methylation could interact with histone modifiers or regulate miRNA to control the growth. Indeed, the importance of miRNA in growth of Nile tilapia has already been shown [68] and some miRNA features (e.g., arm-switching) can be sex-dependent in this species [69].

Several genes related to epigenetic mechanisms of gene expression regulation were also associated with DM CpGs found in males. For example, hdac9 was linked with one intronic CpG hypomethylated in large males, xpo5 related to exonic hypermethylated CpG, and PHD finger protein 20 (phf20 gene) associated with two hypomethylated exonic CpGs. PHF20 is an epigenetic reader involved in histone acetylation [41].

Atg14, a gene found previously related to four CpGs hypomethylated in males when compared to females independently of weight, is also DM between large and small males. Indeed, 20 out of 23 CpGs available in the RRBS dataset (87%) were hypomethylated in large males and situated in the tbp2 (TATA-box binding protein 2) gene body [42] (Figure 8). Either tbp2 or atg14 or both genes could be affected by methylation. Tbp2 is required for transcription initiation of promoters that contain a TATA-box [42]. Thus, if the C-terminal region of Tbp2 is affected by tbp2 methylation, it could have drastic consequences on a large set of TATA-box-containing genes. Although autophagy-related genes have never been studied in Nile tilapia, atg14 should receive particular attention in the future studies related to male-specific growth rate. It suggests an attractive pathway related to self-renewal of cells that can be epigenetically regulated in males but not in females, since the atg14 is highly methylated in the latter regardless of their growth rate.

Uni-parental mouth brooding in cichlids has many physiological and anatomical consequences. Mouth-brooding females have different immune defence and buccal microbiota [70], larger mouth [71], smaller gill size, higher metabolic rate, higher aquatic surface respiration threshold than non-brooding females or males [72]. Indeed, brood protection has an additional energy cost, especially under hypoxia conditions [73]. Mouth brooding time in Nile tilapia is approximately 1–2 weeks, and generally speaking, such a short period of food deprivation should not necessarily induce a strong physiological response in many fish species. However, in the case of Nile tilapia, a longer food restriction occurs exactly at the moment of increased energy demand in females (for example, due to churning behaviour) and has to be considered in addition to the fasting period due to reproduction itself. It is therefore plausible that female Nile tilapia rely on starvation-related autophagy mechanisms during longer periods of fasting, which could be related to the differential CpG methylation of atg14 that we observed between males and females.

Five strongly hypomethylated CpG sites in large males (42% of all available gene-specific CpG sites) were found in the promoter region of the gene encoding for HADH, a marker of fatty acid beta oxidation (catabolic process) occurring in muscle after exercise [74]. Finally, a group of three related genes, gadd45g, gdf3, and ppargc1a, was associated with hypomethylated CpG sites in large males. GADD45G is encoded by a member of a family of genes affected by agents damaging DNA, involved in DNA demethylation [75] or causing growth arrest and muscle atrophy due to stress [43]. Indeed, gadd45 gamma isoform induces the transcription of p21 through the activation of the p38 MAPK (mitogen-activated protein kinase) pathway [76]. An alpha isoform of gadd45 gene has been studied in hybrid tilapia species, although not in muscle tissue and only in the context of bacterial infection [77]. Environmentally induced methylation patterns of gadd45g have been studied in human cells [78] and found to be associated with different carcinomas [79]. Interestingly, the gadd45a gene promoter is transcribed by pparγ gene (peroxisome proliferator-activated receptor gamma) [80,81], regulating gdf3 [82] and closely related to ppargc1a (PGC-1alpha), both hypomethylated in large males. Ppargc1a is associated with diabetes [83], and is considered as a key regulator of energy metabolism [45] and mitochondrial biogenesis [84]. Its decreased expression in muscle is related to abnormal insulin function [85] and hypermethylation of its non-CpG sites is associated with silenced DNA methyltransferase 3B (dnmt3b) and with reduced mitochondrial density [86]. Gdf3 is involved in muscle metabolism [44]. Taken together, these data suggest that epigenetic modulation associated with growth is sex-specific in Nile tilapia. The analysis of the common locations of weight-related DM CpG sites found separately in males and females shows that the total number of sites is not only very low (57), but also that half of them display opposite patterns of growth-specific methylation in males compared in females.

The transcription factor gata6 gene, involved in regulation of smooth muscle growth [47], was associated with exonic DM CpG sites that consistently showed the same patterns of methylation in large fish but was also detected in both growth- and sex-specific comparisons (Figure 7), indicating that gata6 may be involved in transcriptional regulation related to both growth and sex in Nile tilapia.

The importance of sex in growth-related methylation patterns in Nile tilapia is highlighted again by a very low number of DM CpG identified in the sex-independent comparison between large and small fish. Indeed, when large females and large males are grouped together and compared to small fish, only 117 DM CpG sites are found, 10 times less than in the weight-specific comparison between large and small males. Only a few muscle growth-related genes associated with DM CpG sites were unique to the comparison between large versus small fish regardless of sex, such as cdh2 (cadherin 2) [46] (Figure 7).

Conclusions

This study provides a genome-wide reduced representation profiling of cytosine methylation at a single nucleotide level in Nile tilapia muscle. The growth-related methylome in fast muscle differed significantly between males and females, thus confirming our initial hypothesis of sex-related differences in methylation patterns between large and small fish. For example, in males, genes related to stress-induced inhibition of growth rate, fatty acid metabolism, and autophagy were associated with many cytosines showing differential methylation between large and small individuals. However, none of these genes was DM between large and small females, suggesting a sex-specific epigenetic regulation of different gene networks. The autophagy-related gene, atg14, is particularly interesting, since its epigenetic regulation or alternative splicing may affect growth in a sex-dependent manner and whose methylation could be affected during domestication via selective breeding. Several growth-related genes involved in epigenetic regulation of gene expression such as DNA methylation, miRNA, or histone modifications were DM, which indicates that the growth phenotype in Nile tilapia may depend on the interplay between epigenetic mechanisms at different levels involving a large array of biological processes that could be affected during the process of domestication-related growth improvement. In addition to providing novel insight into the epigenetic regulation of muscle growth in fish, the potential epigenetic marks associated with growth rate may be used in marker-assisted selection, thus improving selective breeding and increasing sustainability of fish farming in the future.

Material and methods

Sampling

Nile tilapia eggs were collected from the mouth of wild females caught in the Nile river, Egypt (location GPS: 25°39ʹ56″ N, 32°37ʹ07″ E) and transported to our research station at Nord University (Norway) where they hatched and were reared in a recirculating aquaculture system (pH = 7.6, temperature = 28°C, 11 h dark/13 h light cycle) during 5 months at a density of 27 fish/m³, which corresponded to a maximum density of 6 kg/m³ at the end of the experiment. The fish were fed ad libitum with 0.15–0.8 mm Amber Neptun pellets (Skretting, Sjøhagen 3, 4016 Stavanger, Norway) and kept in the same tank, so as to avoid any potential bias in methylation patterns due to environmental differences (tank effect). They were euthanatized with clove oil (Sigma Aldrich, 3050 Spruce St, St. Louis Missouri, 63103, USA) using a 1:10 mix of 15 mL clove oil with 95% ethanol diluted in 10 L of water and then fast muscle was carefully excised from 24 fish (12 biggest and smallest individuals from each sex, Figure 1) just prior to snap-freezing and then stored at −80°C until DNA extraction. DNA was extracted with the DNEasy Blood & Tissue kit (Qiagen GmbH, Qiagen Strasse 1, 40724 Hilden, Germany) according to the manufacturer’s recommendations. DNA quality check was performed by NanoDrop and DNA quantification was performed with Qubit (Thermofisher Scientific, 168 Third Avenue, Waltham, MA USA 02451) and Tape Station (AgilentTechnologies, 5301 Stevens Creek Blvd, Santa Clara, CA 95051, USA) Assays.

Library preparation

Library preparation for RRBS analysis was performed with the Premium Reduced Representation Bisulfite Sequencing (RRBS) Kit (Diagenode, Rue Bois Saint-Jean, 34102 Seraing (Ougrée), Belgium), essentially following the manufacturer’s instructions. The protocol was optimized by including a double digestion (Taqα1 + MspI) instead of single MspI digestion, which produced more fragments and of more appropriate length (Supplementary Figure 1). In short, normalized DNA samples (100 ng, 26 μL) were sequentially digested with MspI (16 h at 37°C) in CutSmart® Buffer and Taqα1 (8 h at 65°C, followed by deactivation 20 min at 80°C). The digested fragments were then end-repaired, followed by addition of dA. Next, 5ʹ and 3ʹ adapters were ligated, and each sample was size-selected (discarding all reads below 140 bp) with AMPure XP beads (Beckman Coulter, Miami Campus 11800, S.W. 147th Avenue, Miami, Florida 33196, USA), PCR-amplified (13 cycles) and purified again (AMPure XP beads). Finally, the libraries were pooled at 1.6 nM, diluted and single-end sequenced with the NextSeq 500/550 High Output v2 Kit, 75 cycles (Illumina, 5200 Illumina Way, San Diego, CA 92122, USA), according to the manufacturer’s recommendations. Raw reads were converted to FASTQ format, demultiplexed using bcl2fastq V 2.16 (Illumina) and their quality was examined with the FastQC tool (Babraham Bioinformatics). All reads are available at the SRA database under the reference PRJNA506143.

Bioinformatics and statistical analyses

The significance of weight differences between small and large males and females was determined by Student’s t-test, after verifying normality and equal variance assumptions.

Raw sequencing reads were adapter-trimmed with Cutadapt (Babraham Bioinformatics) [87] and aligned to the reference genome ASM185804v2 (orenil3, NCBI, https://www.ncbi.nlm.nih.gov/) using the Bismark v0.19.1 pipeline (Babraham Bioinformatics) [88]. The two strands of oreNil3 were modified in silico (conversion of all C’s to T’s) and indexed according to Bowtie2 [89] requirements. Reads were mapped to the original and in silico-modified genome using Bismark (with parameters: -q – p – N1 – bowtie2). The resulting CpG coverage files were used as input in MethylKit package [90] in R to calculate the coverage, methylation levels of each sample, and differences in methylation among samples. The statistical method used to detect DM cytosines between groups was logistic regression, and the significance thresholds were FDR <0.01 and a minimum difference in methylation level of 25%. Small fish and females were used as reference in size and sex comparisons, respectively. For each CpG site, if the group compared to the reference has a methylation level at least 25% higher, this CpG site was considered as hypermethylated. Conversely, if the methylation level is at least 25% lower in the target than in the reference group, the CpG site was considered as hypomethylated. Integrated Genome Viewer [91] was used to visualize the methylation profile of all CpG from each sample and Circos [92] was used to represent and locate all DM CpG sites. The O. niloticus annotation release 103 was used as input in Homer [93] to locate CpG sites within genomic context (introns, exons, promoters, intergenic regions) and to assign the gene reference in RefSeq format with the closest transcription start site to each DM cytosine. Gene references were then converted to gene symbols and full gene names using the NCBI browser and manually associated with function using the GeneCards browser [94] and the literature. In addition, Functional Enrichment Analysis was performed with DAVID [95] using a custom background gene list corresponding to all the CpG sites analysed in RRBS dataset, and specific to each comparison.

Authors’ contributions

TP carried out the wet lab and bioinformatics analysis, interpreted the results, and wrote the article. SB contributed significantly to wet lab analysis (DNA extraction and RRBS protocol) and revised the article. IK contributed significantly to sampling, to bioinformatics analysis and revised the article. JMOF conceived the study, designed the experiment, provided reagents, and contributed significantly to wet lab analysis, bioinformatics analysis, interpretation of results, and article revision. All authors read and approved the manuscript.

Ethics statement

All procedures involving animals were performed according to the guidelines of the Norwegian Animal Research Authority (FOTS ID 1042) and approved by the Nord University (Norway) ethics committee.

Acknowledgments

We are thankful to Hilde Ribe, Steinar Johnsen, Øivind Torslett and particularly to Kaspar Klaudiussen (Nord University, Norway) for their assistance with fish husbandry.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental material

Supplemental data for this article can be accessed here.

Additional information

Funding

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme [grant agreement no 683210] and from the Research Council of Norway under the Toppforsk programme [grant agreement no 250548/F20].