Early fish domestication affects methylation of key genes involved in the rapid onset of the farmed phenotype.

Animal domestication is a process of environmental modulation and artificial selection leading to permanent phenotypic modifications. Recent studies showed that phenotypic changes occur very early in domestication, i.e., within the first generation in captivity, which raises the hypothesis that epigenetic mechanisms may play a critical role on the early onset of the domestic phenotype. In this context, we applied reduced representation bisulphite sequencing to compare methylation profiles between wild Nile tilapia females and their offspring reared under farmed conditions. Approximately 700 differentially methylated CpG sites were found, many of them associated not only with genes involved in muscle growth, immunity, autophagy and diet response but also related to epigenetic mechanisms, such as RNA methylation and histone modifications. This bottom-up approach showed that the phenotypic traits often related to domestic animals (e.g., higher growth rate and different immune status) may be regulated epigenetically and prior to artificial selection on gene sequences. Moreover, it revealed the importance of diet in this process, as reflected by differential methylation patterns in genes critical to fat metabolism. Finally, our study highlighted that the TGF-β1 signalling pathway may regulate and be regulated by several differentially methylated CpG-associated genes. This could be an important and multifunctional component in promoting adaptation of fish to a domestic environment while modulating growth and immunity-related traits.

Introduction

Animal domestication is a process of modifying a population phenotype through anthropic selection [1] and concerns species chosen for their high value from a human perspective. It starts with wild animals that are kept under controlled conditions and selectively bred or cross-bred for traits of interest, such as improved growth. The environmental conditions during domestication are very distinct of those to which the organism is exposed in the wild. For example, farm-like conditions often include captivity and high density rearing, regular and ad libitum feeding, highly homogenous diet and stable abiotic conditions (temperature, photoperiod), assisted health care and the absence of predatory pressure [2], but also continuous contact with humans. Such drastic changes in the environment can lead to strong phenotypic reaction at an early stage, which, if genetically encoded and repeated over many generations, can eventually give grounds for adaptation. Environmentally-induced changes in phenotype over generations are accompanied and often fostered by strong artificial selection, especially when phenotypic traits sensitive to the environment are also simultaneously targeted by selective breeding [3].

Interestingly, since the farm environments are alike and target traits that are often the same among highly unrelated species, their domestication leads to a set of convergent phenotypic traits, often referred to as the domesticated phenotype [4]. Among the most known traits of domesticated phenotype, we can cite altered behaviours such as increased tameness or reduced stress responsiveness, body proportion and size modifications with reduced locomotory capacities and physiological changes related to increased growth potential [5]. Domestication can also lead to unintentional phenotypic traits, such as decreased immunological resistance [6], a consequence of rearing in relatively sterile conditions or co-selection with targeted traits such a growth [7].

The presence of unwelcome traits in artificially generated domestic phenotype suggests that the process of domestication is not yet fully controlled, and its mechanisms not completely understood. For example, various studies highlight an extremely fast pace of stable phenotypic modifications in a population under ongoing domestication, i.e., within an evolutionary time too short to be explained only by classical genetics. Standing genetic variation can be low in a limited-size founder population F0 [8] and natural occurrence of mutations and subsequent artificial selection can be accounted only in trans-generational intervals, while adaptive changes can occur even within the first generation of organisms submitted to a domestic environment [9,10]. Rapid changes in gene expression can be due to epigenetic mechanisms such as DNA methylation, that not only modulate the response to the new environment but can also be stable over generations. In comparison with genetic variation, epigenetic variation is more likely to have higher rates of spontaneous mutation [11] and to mediate a more sensitive reaction to environmental changes [12,13]. Therefore, in case of environmental changes too fast for an appropriate genetic response, epigenetic mechanisms can have a crucial role in modulating genetically driven phenotypic plasticity and by rapidly generating phenotypic variants that match new environmental conditions and that could be subsequently used in the process of adaptation or artificial selection.

A large spectrum of phenotypic traits associated with epigenetic patterns in domestic animals has been studied in several highly valuable species. For example, DNA methylation has been related to disease resistance [6,14], growth [15], response to diet [16] or rearing conditions [17,18]. These studies raise the hypothesis that this epigenetic mechanism has a rather ubiquitous role in modulating environmentally induced phenotypic traits in farm animals and, in case of commercially valuable traits, could provide a substantial help when establishing selective breeding schemes.

There has been an increasing interest in epigenetics research applied to aquaculture [19], since it has become the fastest growing food production sector in the world [20]. While some fish species such as carp have a known long history of domestication, many new aquatic species are currently at early stages of domestication [21]. Thus, a deeper knowledge of epigenetic mechanisms related to the process of domestication can have a substantial impact on a further development of aquaculture field. In our study, we focused on one of the most important species in aquaculture, Nile tilapia (Oreochromis niloticus), to investigate changes that occur in the fast muscle methylome within a single generation domestication. We expected to detect epigenetic differences underlying fast adaptation to aquaculture conditions and to get an insight into methylation patterns contributing to the onset of ‘domestic phenotype,’ i.e., potential epigenetic marks of fish domestication.

Results

Experimental groups and library characterization

The number of eggs obtained from each female was 192.8 ± 66.1 (mean ± S.D., n = 6) and the overall mortality until the late larval stage was 9.0% ± 2.8% (mean ± S.D., n = 6) (Table S1). After a single generation of domestication, Nile tilapia females from the domestic group were 2-fold heavier than their wild counterparts, despite being of similar age and reproductive status (Figure 1). For the sake of simplicity, the term ‘domestic’ is henceforth used to designate the group of fish undergoing domestication.

Figure 1.

Average weight (g) of 12 females from the wild and domestic groups. Domestic females were statistically heavier than small females (paired t-test, p-value <0.01).

A mean of 25 M raw reads was obtained per library, of which 11 M were uniquely mapped (i.e., 44% mapping efficiency) and an additional 10 M reads could be multiple mapped (Figure 2, Table S2). Neither individual samples nor compared groups differed 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 was around 200-fold (Figure S1). The absence of a peak at the right side of the histogram indicated that the data did not contain PCR duplicates.

Figure 2.

Number of total raw, quality-trimmed, adapter trimmed, uniquely mapped and multiply mapped reads in the domesticated- and wild female group.

Characterization of the fast muscle methylome in Nile tilapia

The global cytosine count was 177 million with a fold coverage of 2.2–2.4 (Table 1). Average methylation corresponded to 60% and did not differ within or between the two analysed groups (Table S3). The mean methylation levels per sample and per chromosome were homogenous, with the exception of chromosome LG3 showing 73% methylation (Table S4). However, the genomic context of the Reduced Representation Bisulphite Sequencing (RRBS) dataset differed from that in the whole Nile tilapia genome (Figure 3). Indeed, intergenic-annotated cytosines were enriched almost two-fold in the reduced representation genome compared to whole genome (Pearson’s Chi-squared, p-value < 2.2e-16). A two-fold decrease in CpG sites within intronic regions was found in the RRBS dataset when compared to the whole genome (Figure 3) (Pearson’s Chi-squared, p-value < 2.2e-16). Different enrichments between reduced- and whole genome cytosines were also observed in exons and TTS (Pearson’s Chi-squared, p-value < 2.2e-16).

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.

Group	Cytosine count	Cytosine fold coverage	%mCG	%mCHG	%mCHH
Domestic	187 M ±1292	2,4	60	0,6	0,4
Wild	167 M ± 857	2,2	61	0,6	0,35

Figure 3.

Genomic context of CpG sites found in the whole genome, our RRBS dataset and associated with hypo- and hypermethylated sites.

Unequal distribution functional relevance of differentially methylated (DM) CpG sites between wild fish and their progeny undergoing domestication

A total of 538,494 CpG positions that were covered at least 10 times and present in every sample at least once were used for analysis of differential methylation between the groups of wild (reference) and domestic group. Interestingly, the genomic context was similar between the sets of hyper- and hypomethylated CpG sites (Figure 3), except for the promoter region, associated in the domestic group with 5 times more hypomethylated sites than the hypermethylated ones (Pearson’s Chi-squared, p-value = 0.002) (Figure 3). Overall, the genomic context of DM CpG sites was similar to the one of all CpGs analysed in this study. In addition, a total of 6 differentially methylated regions (DMR) of 1000 bp (4 hypermethylated and 2 hypomethylated) have been found (Table S5), containing a minimum of 10 CpG each and each covered at least 3 times.

Out of a total number of 715 differentially methylated CpG sites (FDR < 0.01) (Figure 4, Table S6), the number of hypo- and hypermethylated CpG sites was similar and represented on average 0.7‰ of all CpG sites analysed. The chromosomal distribution differed between hypo- and hypermethylated CpG sites (Table S7). For example, chromosome LG3 was represented by a similar number of hypo- and hypermethylated CpG sites, but chromosome LG18 contained 3 times more hypomethylated sites than hypermethylated ones. Within-chromosomes, the distribution of DM CpG sites did not show any specific pattern, except for the longest chromosome LG23, which contained all the DM CpG sites within the first 50% of its sequence (Figure 5). The lack of DMCpG sites in the second half of chromosome LG23 can be explained by a low RRBS coverage of this region (only 1664 CpGs analysed were located further than 44 M nucleotides) compared to the first half of LG23 (28,310 CpGs located within the first 44 M nucleotides of LG23). The Principal Component Analysis on all cytosine levels available in the RRBS dataset showed a clear separation between wild and domestic groups (Figure 6, Table S8).

Figure 4.

Volcano plot representing q-value and methylation difference in all CpG sites compared between wild and domesticated fish. Red dots represent differentially methylated CpG sites (methylation difference >25%, FDR <0.01). Several genes associated with most extreme CpG methylation differences are indicated. In addition, blue dots highlight all the 15 DM CpG sites associated specifically with stub-1 gene.

Figure 5.

Chromosomal distribution of differentially methylated CpG sites between wild females and their progeny undergoing domestication (716 positions, FDR < 0.01, methylation difference >25%). Wild females are set as reference. Histograms pointing inwards (dark green) and outwards (dark red) represent hypomethylated and hypermethylated sites in domestic group, respectively.

Figure 6.

Principal component analysis of 12 fish separated in the domestic (d, blue) and wild group (w, red).

The functional annotation only showed a significant enrichment of functions related mostly to rRNA metabolism and modifications (Table S9). These limited results are due, at least in part, to the relatively small number of genes associated with DM CpGs. Among the 489 manually annotated genes associated with a total of 715 differentially methylated CpG sites (Table S6), we focused on genes linked to metabolism, muscle growth, immunity and epigenetics, since they are directly relevant to the domestication phenotype (Table 2).

Table 2.

List of functionally relevant selected genes associated with DM CpG sites in Nile tilapia females. Six domesticated individuals were compared to six wild individuals. Only the results with q-value < 0.01 and methylation difference >25% are retained. Positive percentage values indicate hypermethylation in domesticated individuals when compared to wild 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.

Chr	Start	Diff. (%)	FDR	Strand	Annotation	Distance Nearest TSS	RefSeq Nearest Promoter	Gene Symbol	Function	DM CpG count	Frequency Ratio
chrLG22	5,571,043	−31,95	2,7E-12	-	intron	28,408	XM_005457545.4	clcn1	muscle functioning	2	0,17
chrLG22	5,570,993	−30,88	3,94E-10	-		28,458
chrLG13	2,252,388	26,39	1,07E-09	-	intergenic	40,827	XM_005462223.4	cst7	immunity	13	0,10
chrLG13	2,252,441	25,26	1,32E-14	-		40,774
chrLG13	2,252,449	29,88	1,39E-19	-		40,766
chrLG13	2,252,466	26,09	2,03E-15	-		40,749
chrLG13	2,263,034	−25,32	2,56E-13	-		30,181
chrLG13	2,252,456	30,32	3,24E-19	-		40,759
chrLG13	2,236,611	31,90	3,53E-19	-		56,604
chrLG13	2,263,001	−28,05	3,76E-07	-		30,214
chrLG13	2,263,005	−26,22	3,79E-14	-		30,210
chrLG13	2,252,461	30,22	5,34E-21	-		40,754
chrLG13	2,263,013	−28,09	5,37E-17	-		30,202
chrLG13	2,263,003	−29,85	7,27E-13	-		30,212
chrLG13	2,252,374	31,93	8,92E-13	-		40,841
chrLG11	24,202,430	36,90	1,02E-38	+	exon	−15,724	XM_025911682.1	cic		1	0,01
chrLG4	33,455,755	26,60	1,25E-04	+		−7504	XM_019357563.2	stub1	muscle growth/ immunity	15	0,25
chrLG4	33,455,772	25,26	3,57E-04	+		−7487
chrLG4	33,455,758	25,18	3,70E-04	+		−7501
chrLG4	33,456,304	29,82	1,10E-15	+		−6955
chrLG4	33,456,328	29,87	1,13E-15	+		−6931
chrLG4	33,456,343	25,67	1,57E-11	+		−6916
chrLG4	33,456,336	27,73	1,65E-13	+		−6923
chrLG4	33,456,356	27,75	1,66E-13	+		−6903
chrLG4	33,456,299	27,63	1,88E-13	+		−6960
chrLG4	33,455,740	29,05	1,97E-05	+		−7519
chrLG4	33,456,318	28,54	2,47E-14	+		−6941
chrLG4	33,456,365	29,30	4,09E-15	+		−6894
chrLG4	33,456,339	30,21	4,19E-16	+		−6920
chrLG4	33,456,368	26,41	5,85E-05	+		−6891
chrLG4	33,455,761	27,13	9,37E-05	+		−7498
chrLG23	18,304,157	−26,81	1,37E-12	-	intergenic	−5091	XM_025903147.1	acvr2a		3	0,06
chrLG23	18,304,202	−26,58	3,32E-12	-		−5136
chrLG23	18,304,207	−26,33	5,60E-12	-		−5141
chrLG11	38,128,296	25,62	2,17E-12	-		−22,285	XM_025911147.1	smad4		2	0,03
chrLG11	38,128,308	32,40	6,16E-16	-		−22,297
chrLG5	32,577,336	25,41	1,01E-18	-	exon	18,307	XM_005450144.4	pdzrn3		1	0,01
chrLG7	144,049	28,36	1,94E-08	+		533	XM_013266159.3	dgkz		9	0,06
chrLG7	144,116	27,89	3,99E-05	+		600
chrLG7	144,085	25,96	4,50E-07	+		569
chrLG7	144,064	25,80	5,61E-07	+		548
chrLG7	144,107	25,68	6,67E-07	+		591
chrLG7	144,056	25,51	7,83E-07	+		540
chrLG7	144,077	25,54	8,27E-07	+		561
chrLG7	144,104	25,51	8,82E-07	+		588
chrLG7	144,061	25,43	9,18E-07	+		545
chrLG15	1,967,239	32,38	1,63E-15	-	intergenic	−2028	XM_003440691.4	nfkbia		1	0,03
chrLG13	33,423,455	26,51	1,52E-25	+	intron	−53,177	XM_025897472.1	itbp1		1	0,02
chrLG18	33,712,851	−29,45	1,29E-07	-	exon	6984	XM_005457639.2	dusp28	muscle growth	2	0,22
chrLG18	33,712,834	−27,79	6,98E-07	-		7001
chrLG14	37,954,443	31,27	4,06E-14	-	intergenic	−31,600	XM_013273630.2	shkbp1		1	0,01
chrLG7	56,727,304	−25,94	2,16E-14	+		30,107	XM_005470918.4	sirt3		1	0,02
chrLG18	34,613,693	−28,91	1,30E-08	+	exon	666	XM_025900193.1	cdyl	muscle growth/ epigenetic	2	0,07
chrLG18	34,613,695	−28,00	1,78E-08	+		668
chrLG15	26,870,849	−28,34	2,90E-09	-	intergenic	41,311	XM_025898869.1	ubr7	epigenetic	1	0,01
chrLG9	13,162,946	−25,33	1,39E-09	+	exon	26,977	XM_003439271.5	fbxl7	autophagy/ epigenetic	1	0,06
chrLG6	12,633,677	26,34	9,52E-09	-	intron	−24,630	XM_003455961.5	alkbh5		1	0,02
chrLG20	24,778,943	−29,90	1,12E-12	-	intergenic	−39,699	XM_005478335.4	dcp1a		1	0,06
chrLG2	20,651,068	−25,62	1,01E-09	-	exon	9293	XM_005467513.4	atg2a	autophagy	1	0,20
chrLG19	26,567,932	28,38	4,18E-12	+	intergenic	25,554	XM_005477597.4	capn2		1	0,02
chrLG16	32,968,330	27,11	5,79E-06	+		29,135	XM_025903828.1	itb2		2	0,13
chrLG16	32,968,309	34,14	1,32E-10	+		29,114
chrLG3	24,036,571	−31,27	7,28E-27	+		−8592	XM_003447550.4	capn1		1	0,05
chrLG18	27,326,614	26,18	4,90E-05	-		−20,496	XM_019347353.1	acot1	obesity/ high fat diet	1	0,02
chrLG13	2,031,574	47,23	9,08E-40	+		−39,871	XM_019366617.2	elovl6		1	0,01
chrLG1	30,667,269	27,02	2,82E-03	+	exon	1920	XM_005447642.4	irx3		3	0,11
chrLG1	30,667,250	44,49	1,94E-08	+		1901
chrLG1	30,667,254	32,75	9,12E-05	+		1905
chrLG11	1,390,538	−30,32	9,72E-06	-	intron	−5094	XM_005459395.4	adcy3		1	0,03
chrLG7	1,438,772	25,81	1,96E-09	-	intergenic	26,265	XM_019360931.2	mlkl		1	0,04
chrLG10	32,876,536	26,91	4.52e-05	-	exon	1434	XM_003457007.5	egr1	muscle growth obesity/ high fat diet	1	0,01

DM CpGs associated with muscle-growth related genes

Two hypomethylated exonic CpGs were associated with dusp28, which encodes a protein from a family of dual-specificity phosphatase (DUSP) that can activate the p38 mitogen-activated protein kinase (MAPK) signalling pathway [22] and therefore regulate muscle stem cells [23]. Three hypomethylated DM CpG sites were linked to acvr2a, encoding Activin receptor type 2a that is linked to myostatin activity in muscle growth [24]. Another hypomethylated CpG site was associated with sirt3, a gene involved in metabolic flexibility of skeletal muscle (i.e., ability to switch between glucose and lipid oxidation) by impacting downstream signalling through peroxisome proliferator–activated receptor (PPAR)γ coactivator-1α [25]. The activity of sirt3 is regulated by nutrient signals and contractile activity of the muscle [26] and has been linked to the fast growing fish of gilthead sea bream species [27]. Hypermethylated CpG sites in the domestic group were present in pdzrn3, encoding a ubiquitin ligase regulating myoblast differentiation [28] and in shkbp1, encoding a SHK3BP1 binding protein involved in development of skeletal muscle fibres [29]. Moreover, hypermethylation was also found in an exonic CpG site related to egr1, an early growth response gene. This transcription factor regulates myoG gene expression and promotes differentiation of muscle satellite cells [30], but is also affected by high intensity exercise [31]. Egr1 also inhibits another transcription factor, MEF2, important in muscle differentiation [32]. Importantly, gene body methylation of egr1 was associated with intramuscular fat in pigs at different developmental levels [33], its expression is up-regulated by myoD [34] and downregulated by Transforming growth factor-β1 TGF-β1 [35]. Two hypomethylated intronic CpG sites were also associated with clcn1, coding for Chloride channel protein 1, which, although not involved in muscle growth specifically, plays a role in muscle physiology through membrane repolarization after muscular contraction [36].

DM CpG associated with both muscle growth and immunity-related genes

Several DM CpG associated genes had relevant functions not only in muscle growth but also in the immune response (Table 2). For example, stub1 was associated with the highest number of DM CpGs, all hypermethylated in domestic fish and situated in exons. There was a clear difference in methylation levels between wild and domestic fish in individual CpG sites associated with stub1 (Figure S2), which are present in two close clusters in a volcano plot (Figure 4). Stub1-associated DM CpG were also detected by DMR analysis (Table S5). This gene codes for E3 Ubiquitin ligase CHIP (carboxyl terminus of Hsc70-interacting protein), which is involved in ubiquitin-mediated protein degradation and regulation of TGF-β1 [37]. TGF-β kinase has a well-known function in inflammation and wound healing [38], and is important for muscle regeneration after damage or exercise [39]. TGF-β also induces skeletal muscle atrophy [40] by delaying myoblast differentiation while increasing cellular proliferation [41]. Another gene associated with hypermethylated CpG was ltbp1, encoding the Latent transforming growth factor beta binding protein 1, which is critical for the assembly of TGF- β1 [42], a cytokine important in muscle regeneration [43]. Interestingly, a significant correlation between gene expression of ltbp1 and the methylation level of its CpG sites was found in sheep [44]. Another gene known as a key regulator of TGF-β1 signalling pathway and associated hypermethylated CpG sites was smad4. This gene is also affected by diet in grass carp and is associated with muscle hardness through regulation of type 1 collagen expression [45]. Smad4 is also involved in skeletal muscle regeneration by promoting the expansion of satellite cell derived progenitors [46]. Besides the muscle growth-related function, this transcription factor regulates the activity of PPARγ via TGF-β1 [47], but also plays a direct role in the immune response [48]. Nine hypermethylated exonic CpG sites were related to dgkz encoding diacylglycerol kinase zeta (DGKZ), which is involved in cardiomyocyte hypertrophy through the PPAR (peroxisome proliferator-activated receptor)-DGK pathway [49]. This gene was associated with feed conversion ratio in broilers by GWAS study [50] and its methylation was suggested as a biomarker of prostate cancer [51]. Importantly, roles of dgkz in myogenic differentiation [52], in muscle hypertrophy and myoblast fusion through mTOR-dependent signalling [53,54] and in T cell activation [55] were experimentally validated, and multiple studies of its function in modulating lipid metabolism has been reviewed [56]. Two hypermethylated exonic CpGs were associated with the gene itb2 encoding Integrin beta-2, with known functions in skeletal muscle hypertrophy [57] and immune response [58]. Another DM CpG-related gene was nfkbia, which codes for an inhibitor of the transcription factor NF-kappa-B, the latter being important in muscle atrophy [59] and inflammation [60].

Immunity-related genes associated with DM CpG

Thirteen hypo- and hypermethylated exonic CpG sites in domestic group were associated with cst7 encoding Cystatin F, a cysteine protease inhibitor involved in eosinophil survival [61] and in immune response [62]. Another gene associated with a hypermethylated exonic CpG was the gene cic encoding Capicua, a transcriptional repressor that maintains peripheral immune tolerance via T-cell homoeostasis [63] and whose methylation is inversely correlated with its transcription level in pigs [64].

Genes associated with DM CpG and related to epigenetic mechanisms

Several DM CpG sites were associated with genes that have important roles in histone binding, RNA degradation and RNA methylation. For example, ubr7 encoding E3 Ubiquitin ligase UBR7, is associated with histone marked segments of the genome [65], interacts with histone H3 [66] and has been suggested to function as histone E3 monoubiquitin ligase [67]. Two hypomethylated exonic CpGs were associated with the gene cdyl encoding the chromodomain protein CDYL, which can bind to H3K27 histone and interact with Polycomb Repressive Complex 2 [68], a histone methyltransferase that is also involved in regulation of muscle cell differentiation [69]. Another DM CpG-related gene was dcp1a coding for an mRNA decapping enzyme involved in RNA degradation and translation regulation [70]. Interestingly, the closely related gene dcp2 is also involved in epigenetic regulation of autophagy [71]. Finally, a gene alkbh5 associated with hypermethylated intronic CpG encodes for a m6A RNA demethylase with epigenetic activity linked to autophagy [72] and to the transcription factor EB, involved in metabolic flexibility [73]. Moreover, alkbh5 expression can be affected by diet supplementation [74]. These last two genes were not only related to epigenetic processes but also to autophagy.

Genes associated with DM CpG and related to autophagy and proteolysis

Hypomethylated CpG sites in domestic fish were related to fbx7, atg2a and capn2. Fbxl7 encodes for F-box and leucine rich repeat protein 7 that modulates mitochondrial functions such as proteasomal degradation by negatively regulating an anti-apoptotic protein survivin [75]. Atg2a, autophagy related gene 2, is involved in glucose starvation related autophagy with a suggested function in controlling the extent of autophagosome membrane formation [76,77]. Capn2, encoding Calpain 2, is involved in autophagosome formation and protein degradation [78] and was associated with meat quality in birds [79]. Interestingly, the expression of this protease is positively regulated by testosterone in fish [80]. One hypomethylated CpG was related to the gene capn1 encoding Calpain 1, which plays a crucial role in autophagy via autophagy-related genes [81] and is associated with nutritional state in halibut [82] and gilthead sea bream [83]. The activity of this protease was also induced by a moderate and sustained exercise in the latter [84].

Genes associated with DM CpG and related to obesity and high fat diet

Several hypermethylated CpG sites were associated with genes related to obesity and fat composition, such as adcy3, mlkl, irx3, elovl6 and acot1. The expression of adcy3 increases in skeletal muscle in response to fat overfeeding [85] and its related protein is associated with diet-induced obesity [86]. Mlkl codes for mixed lineage kinase domain-like protein, also involved in obesity-induced metabolic complications [87]. Three exonic DM CpG sites were associated with irx3, encoding for Iroquois homeobox 3, strongly related to obese phenotype and obesity-associated gene fto [88]. The gene elovl6 encodes a Fatty acid elongase and is associated with a QTL effect on fatty composition [89] and obesity-induced insulin resistance in vertebrates [90]. Acot1 codes for Acyl-coenzyme A thioesterase 1, which is upregulated in response to high fat overload in animals fed with high-fat diet [91].

Discussion

In this study, we report the reduced representation profiling of methylated cytosines in fast muscle between wild females and their hatchery-born and reared female offspring, i.e., subjected to one generation of domestication. The overall mortality values were relatively low and within the range expected for Nile tilapia (Table S1), indicating that the effect of artificial selection was negligible. Therefore, our study gives insight into epigenetic patterns influenced by the early stages of domestication. Little is known on the genetic differences between the fish, however, the experimental design and analysis (same sampling site for the wild fish, parent-offspring relation of each pair of fish used as covariate) allowed to rule out the influence of genetic background in the study, if any. Moreover, a study performed in coho salmon has shown that hatchery conditions have stronger impact on the epigenome than fish origin and that no significant genetic differentiation arose just after one generation of domestication process [92]. A recent study in sea bass under early stages of domestication has shown that domestication-related methylation patterns not only could explain the onset of domestic traits, such as lower jaw malformations, but also that differentially methylated sites correspond to genetic polymorphisms observed after long periods of selective breeding and coincide with genes under positive selection during domestication process [93].

To the best of our knowledge, this is the first study providing a cytosine resolution methylation profiling in Nile tilapia between generations in the context of early domestication. Previous studies have shown cytosine methylation differences in fast muscle associated with Nile tilapia growth and sex in a single generation of hatchery conditions [94], regional cytosine methylation differences in fast muscle associated with sex in a Nile tilapia hybrid [95] or associated with sex in gonads [96]. A limitation of these studies and indeed the present paper is the use of complex tissues composed by multiple cell types. Differences in tissue composition can result in the identification of additional differentially methylated genes between groups, namely immune-related genes due to the presence of blood or neuronal markers due to variations in innervation. Recent advances in single-nucleus methylome sequencing will make it possible to address this issue in the near future.

Non-tilapiine methylome profiling within the context of domestication has been performed in a few other species, namely in fast muscle of coho salmon [92] and in the steelhead fin tissue [97], sperm and red blood cells [18]. We compared our results with the closest tissue- and context-specific studies [92,94] and found a few common genes associated with differentially methylated cytosines or regions. Genes from a family encoding dual specificity phosphatases (DUSP) and protocadherins (CADH) were associated with cytosine hypomethylation of domesticated groups in both our and the previous study [92], and those encoding for serine/threonine phosphatase kinases (SGK) and syntaxins (STX) were associated with cytosine hypomethylation of domesticated fish in our study as opposed to hypermethylation in the above mentioned study. However, it is difficult to interpret since none of the gene families concerned strictly the same gene or species, and the genes were associated with large differentially methylated regions in coho salmon, in contrast to our single cytosine-specific association. More common genes were found when compared to the single cytosine methylome analysis between fast- and slow growing Nile tilapia [94]. For example, within the family of autophagy-related genes, where atg14 was recently shown as associated with the highest number of growth- and sex-related DM CpG sites. It strengthens the idea that autophagy might be an epigenetically regulated and both growth-specific and sex-dependent mechanism in Nile tilapia.

Although the average methylation level of the Nile tilapia genome is similar to that reported in previous studies, the proportion and distribution of DM CpG sites throughout the genome suggests a non-random pattern of methylation associated with domestication. The systematic analysis of functional ontology of genes associated with these patterns has highlighted the GO terms related to ribosomal RNA modifications, and specifically rRNA methylation. Such rRNA modification can have various consequences on the phenotype, since ribosomes are a general and essential element of protein synthesis. DNA methylation patterns involved in rRNA methylation make a link between pre- and post-transcriptional epigenetic marks and add an extra layer of complexity of epigenetic mechanisms involved in domestication. Moreover, although more processes were not detected as significantly enriched by systematic analyses, we found many DM CpG-associated genes related to growth, immunity, autophagy, diet or epigenetic mechanisms (Figure 7) that can help explaining the changes occurring during the early stages of domestication. One of the main differences between any wild and captive population is the access to food, restrained in the first and often ad libitum in the second group. Farmed fish receive more food at regular times and without wasting energy for foraging, and the food is nutritionally optimized to yield high growth rates. Thus, such important dietary changes can be reflected in gene regulation related to diet components, which in our case were associated particularly with high fat nutritional profiles. The growth-promoting farming conditions were reflected by several DM CpG associated genes, many of which, intriguingly, were also implicated in immunity-related processes. This growth-immunity relationship is in line with the life-history theory predicting trade-offs between physiologically costly functions and it has been shown that, in case of growth-selected fish [98] and other farmed animals [99], energy investment in growth is associated with decreased energy allocation in immunity-related functions. It is difficult to predict the consequences of methylation of these genes in our study, however, it is important to highlight that even when energy intake is not restrained by resource availability, the pleiotropic nature of the epigenetically affected genes can lead to important changes in both traits simultaneously. We also do not rule out the hypothesis that in highly controlled farm conditions, the relaxed pathogen pressures lead first to changes in immune profiles of fish, which in consequence will also affect the growth rate. Several selected genes were also associated with proteolysis and autophagy. Protein lysis might be linked to catabolism, i.e., muscle degradation in order to extract energy when food is scarce. In a context of unrestrained food availability together with regular and frequent feeding periods in farm conditions, methylation of these genes could act as a protective mechanism against catabolism, thus promoting once again muscle growth. Moreover, autophagy related genes have been already identified in this species as growth and sex-related in the context of early domestication [94]. Finally, several differentially methylated genes were involved in other epigenetic mechanisms, i.e., mRNA degradation, histone modifications or RNA methylation, suggesting that early domestication involves changes at different levels of gene-phenotype axis.

Figure 7.

Schematic representation of genes that are associated with differentially methylated CpGs (domestic VS wild females) and whose functions are related muscle growth, autophagy, immunity, diet or epigenetic regulation processes.

Several DM CpG associated genes were related to the Transforming growth factor beta 1 (TGF-β1) signalling pathway and regulation (Figure 8). For example, ltbp1 and stub1 are known to promote or inhibit TGF-β1, respectively, while egr1 and smad4 can be repressed or indirectly involved in repression and induction of other genes via TGF-β1. Downstream regulated genes such as mef2, myoG and ppargy, are known for their key role in myogenesis. Even though this kinase was not detected as associated with differentially methylated cytosines, it seems to be affected by upstream regulators and to involve downstream components that are also associated with DM CpGs. We suggest that in the context of early fish domestication, TGF-β1 signalling pathway is one of the candidates that should receive further attention in functional studies.

Figure 8.

Schematic representation of a gene network associated with differentially methylated CpGs (in blue) and whose functions are related to transforming growth factor beta 1 or to other genes relevant to the process of muscle growth (in black). Red/green and continuous/dashed arrows indicate positive/negative and direct/indirect regulation, respectively.

Taken together, our study shows the importance of the epigenetic changes that occur during the early stages of domestication. The differentially methylated genes observed between the wild and first-generation farmed fish could explain, to some extent, the onset of a ‘domestic phenotype’ in aquaculture. Myogenic and immune functions were common to many pleiotropic genes, while other genes are known to be affected by diet and fat in particular. Differential methylation in autophagy- and proteolysis-related genes can potentially decrease catabolism in an environment of regular feeding and unlimited food availability. We suggest that the modulation of phenotypic traits related to growth and immunity might occur already within the first generation of farmed fish, i.e., very early in domestication, and in the absence of artificial selection, this epigenetic mechanism plays a crucial role in early domestication.

Materials and methods

Ethics statement

This study was approved by the Nord University (Bodø, Norway) ethics board and all procedures involving animals were performed according to the instructions of the Norwegian Animal Research Authority (FOTS ID 1042).

Collection of wild specimens, fish husbandry and sampling

Fertilized eggs were collected from six wild mouthbrooding females in the river Nile, Egypt (location GPS: 25°39'56”N, 32°37'07”). They were then euthanized as detailed below and measured. Their fast muscle was sampled, while their fertilized eggs were transported to our research station at Nord University (Norway), where they hatched and were reared in separate tanks in a recirculating aquaculture system (pH = 7.6, temperature = 28°C, 11 h dark/13 h light cycle, density = 27 fish/ m3, food: Amber Neptun Skretting® feed 0.15–0.8 mm) during 5 months. Then one female offspring corresponding to each wild female was euthanized and its fast muscle sampled.

The wild and domestic fish were euthanized and sampled following the same procedure. In short, they were euthanized with clove oil (Sigma Aldrich, USA) using a 1:10 mix of 15 mL clove oil with 95% ethanol diluted in 10 L of water, the fast muscle was carefully excised just prior to snap-freezing and then stored at −80°C until DNA extraction. DNA was extracted with the DNEasy Blood & Tissue kit (Qiagen, Germany) according to the manufacturer’s recommendations. DNA quality check was performed by NanoDrop and Tape Station (Agilent, USA) Assays, while DNA quantification was performed with Qubit (Thermofisher Scientific, USA).

Preparation of reduced representation bisulphite sequencing (RRBS) libraries

Library preparation for RRBS was performed with the NuGen ovation RRBS methyl-seq system 1–16 (Tecan Genomics, Inc, Redwood City, USA) following the producer’s instructions. Genomic DNA was digested with MspI at 37°C for 1 hr, followed by adapter ligation and final repair. For the bisulphite conversion of adapter ligated libraries and subsequent cleaning, the EpiTect fast bisulphite conversion kit (Qiagen, Hilden, Germany) was used following the manufacturer’s protocol. The resulting bisulphite converted libraries were amplified for 12 PCR cycles and later purified using reagents and recommendations of the NuGen RRBS kit. Quality and quantity of the RRBS libraries were assessed using the TapeStation. Single-end 75 bp sequencing was performed for a pool of 12 libraries on an Illumina NextSeq instrument (San Diego, USA) with 4% PhiX DNA (Illumina) as an internal control, following the instructions for RRBS sequencing from NuGen. Sequencing data are publicly available at NCBI under the Bioproject PRJNA661533 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA661533).

Bioinformatics and statistical analyses

Raw reads were adapter-trimmed with the Nugen-provided script and Cutadapt (Babraham Bioinformatics) [100], and aligned to the latest reference genome O_niloticus_UMD_NMBU (https://www.ncbi.nlm.nih.gov/assembly/GCF_001858045.2/), using the Bismark v0.19.1 pipeline (Babraham Bioinformatics) [101]. The two strands of oreNil3 were modified in silico (conversion of all C’s to T’s) and indexed according to Bowtie2 [102] 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 [103] in R to calculate the coverage, methylation levels (i.e., number of methylated counts over a total count at each site) of each sample and differences in methylation among samples. The statistical method used to detect differentially methylated cytosines (DMC) and regions (DMR) between groups was logistic regression (wild group as control, groups as dependent variable, methylation as independent variable), and the significance threshold were FDR < 0.01 and minimum difference in methylation level of 25%. For DMC, only cytosines covered at least 10 times were considered, and for DMR tiling windows were set at 1000bp, initial per base coverage of 3 and a minimum number of 10 cytosines per region. Circos [104] was used to represent and locate all differentially methylated CpG sites. The Oreochromis niloticus annotation release 104 was used as input in Homer [105] to locate CpG sites within their genomic context (introns, exons, promoters, intergenic regions) and to assign the gene reference in RefSeq format with the closest TSS (transcription start site) to each differentially methylated 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 [106] and the literature. In addition, Functional Enrichment Analysis was performed with DAVID [107] and the top 10 pathways were reported, using a custom background gene list corresponding to all the CpG sites analysed in the RRBS dataset, and specific to each comparison.