Fish as Model Systems to Study Epigenetic Drivers in Human Self-Domestication and Neurodevelopmental Cognitive Disorders

Modern humans exhibit phenotypic traits and molecular events shared with other domesticates that are thought to be by-products of selection for reduced aggression. This is the human self-domestication hypothesis. As one of the first types of responses to a novel environment, epigenetic changes may have also facilitated early self-domestication in humans. Here, we argue that fish species, which have been recently domesticated, can provide model systems to study epigenetic drivers in human self-domestication. To test this, we used in silico approaches to compare genes with epigenetic changes in early domesticates of European sea bass with genes exhibiting methylation changes in anatomically modern humans (comparison 1), and neurodevelopmental cognitive disorders considered to exhibit abnormal self-domestication traits, i.e., schizophrenia, Williams syndrome, and autism spectrum disorders (comparison 2). Overlapping genes in comparison 1 were involved in processes like limb morphogenesis and phenotypes like abnormal jaw morphology and hypopigmentation. Overlapping genes in comparison 2 affected paralogue genes involved in processes such as neural crest differentiation and ectoderm differentiation. These findings pave the way for future studies using fish species as models to investigate epigenetic changes as drivers of human self-domestication and as triggers of cognitive disorders.


Introduction
Domestication is defined as an "evolutionary process that substantially reshapes the genetic, physiological and behavioral profile of a species to adapt to a human-made environment" [1]. Historically and contemporarily, this process has affected the evolutionary trajectories of several economically and culturally important vertebrate species. Domestication usually involves selection of less aggressive and more human-tolerant individuals, plus some other specific features of interest such as meat, wool, and milk production. Nonetheless, new additional phenotypic traits emerge repeatedly in independent vertebrate domestication events, even at the early stages of living in a human-made environment prior to deliberate selection; a phenomenon first noticed by Darwin himself and currently known as the "domestication syndrome" [2]. The domestication syndrome has been predominantly described in mammals, likely due to the large number of mammalian domesticates with a long domestication history, sometimes dating back millennia (e.g., dogs). Phenotypic traits of the domestication syndrome, which are not necessarily behavioral domesticated traits in humans with neural crest development and migration processes. Therefore, cognitive disorders and the gene networks associated with them may be used as models for further testing the HSD hypothesis.
Domestication is a process of adaptation to a new selective environment and has been considered likely to involve epigenetic changes [21][22][23][24][25]. Epigenetic mechanisms offer a way for novel phenotypes to emerge rapidly in response to environmental changes and to prime the offspring, when inherited, to face environments based on the parental experience [26][27][28].
In the first stages of domestication, which coincides with the emergence of domestication syndrome traits, epigenetic changes established during early development might regulate gene expression in the neural crest, and be maintained throughout adulthood and inherited by the offspring. Multigenerational epigenetic inheritance is nearly ubiquitous in diverse animal species (see [29] for review). The persistence of the domestication environment, together with the stability and small effect of epigenetic changes in mild developmental deficits of the neural crest, could be expected to accelerate the adaptation [30]. After several generations, epigenetic changes could be genetically assimilated as genetic variants [25,31,32], hardwiring these changes. Partial evidence for this process comes from studies on mammals (dogs-wolves [23]), birds (red jungle fowl-modern chickens [33]), and fish species. The same process could be hypothesized to account for the first steps of HSD, as most differences between extinct hominins and AMHs are epigenetic by nature, having impacted features that are related to the domestication syndrome, particularly the morphology of the face [34]. To date, DNA methylation, but also DNA hydroxymethylation, have been the epigenetic mechanisms of focus, even though other epigenetic mechanisms like histone modifications and variants or non-coding RNAs may also play a role in domestication. Unfortunately, the latter cannot be studied in ancient DNA remains.
Domestication of fish species has a distinct history from terrestrial vertebrates [35], although it is scientifically considered to represent a similar process [36]. Until the 20th century, the majority of seafood relied on wild animal captures, with a few exceptions like the common carp (Cyprinus cyprio) in China~8000 years ago or Nile tilapia (Oreochromis nilocitus) in Egypt~3500 years ago [36,37]. In the last century, the domestication of aquatic species has expanded rapidly, with an estimated number of 368 vertebrates that have been domesticated for aquaculture, including teleost fish, frogs and reptiles [38]. Nevertheless, the majority of species are at the early stages of domestication, without closed life cycles in captivity and in the absence of deliberate selection for specific traits [38]. Nonetheless, in parallel with the domestication process, phenotypic traits involving the domestication syndrome, with changes in growth, reproduction, morphology, pigmentation and behaviour, have become manifested in domesticated fish [39][40][41]. Furthermore, sequencing of fish genomes has revolutionized vertebrate comparative genomics and has greatly contributed to our understanding of selection targets, evolutionary changes and speciation. Subsequently, fish have been suggested to serve as suitable models for human biomedical research [42,43]. Recently, epigenetic patterns emerging during the first stages of domestication, in the absence of genetic differences, have been studied in salmonids [44,45], European sea bass (Dicentrarchus labrax) [39], Nile tilapia (Oreochromis niloticus) [46,47] and grass carp (Ctenopharyngodon idellus) [48]. These epigenetic patterns of domestication are present in the sperm of several species, i.e., salmonids [45,[49][50][51], showing the potential of intergenerational transfer, while in the European sea bass~20% are found in early embryos, showcasing the importance of developmental aspects during early domestication [39]. Taken together, (1) the recent domestication events in fish, (2) the high degree of parallelism between fish and human domestication, particularly, the absence of deliberate selection in both domestication events, (3) the increasing availability of fish methylomes and (4) the use of fish as animal models in biomedical research, render fish promising candidate models to identify the epigenetic mechanisms that lead to the emergence of HSD, including their abnormal manifestation in neurodevelopmental disorders.
Comparative epigenomic studies between domesticated animals and humans are expected to demonstrate parallel or contrasting processes operating in addition to traditional genetic aspects [15]. Here, we argue that fish hold great advantages as models to study epigenetic drivers in HSD. To test our argument, we use comparative approaches to epigenomic patterns, exemplified by the best-studied modification, DNA methylation, between humans and the European sea bass. The European sea bass was chosen because (1) 25 years of selective breeding resulted in selective sweeps in genes similar to those found under positive selection in all domesticates tested, e.g., glutamate receptors [52,53], (2) it presents traits of the domestication syndrome shared with those found in terrestrial vertebrates, e.g., depigmentation and cranial changes [39], and (3) epigenetic patterns of domestication have been assessed in four tissue types representative of all three embryonic layers, thus reducing bias due to tissue-specificity [39]. In the present study, we compare epigenetic patterns of domesticated sea bass with epigenetic patterns of (1) AMH as opposed to archaic hominins (Neanderthals and Denisovans), and (2) neurodevelopmental cognitive disorders with an abnormal presentation of traits parallel to the domestication syndrome (SZ, WS and ASD; Figure S1). The goal of these comparisons was to investigate whether genes or pathways were consistently altered during the steps of early domestication in European sea bass and humans, with a potential impact on our species-specific distinctive cognition and behavior.

Data Collection
Comparative epigenomic analyses were divided into two major groups including early domesticates of the European sea bass vs. (1) AMH and (2) neurodevelopmental cognitive disorders. For this, we compiled five lists of genes identified as differentially methylated in the literature ( Figure S1).

European Sea Bass Early Domesticates
In European sea bass, we previously conducted work to generate genome-wide DNA methylation patterns (Reduced Representation Bisulfite Sequencing, RRBS) in fish captured in the wild vs. offspring of wild fish reared in hatchery [39]. DNA methylation data from the brain, muscle, liver and testis can be accessed through the NCBI Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/, accessed on 30 May 2022) with accession codes GSE104366 and GSE125124. Since these data were published, the European sea bass genome has been included in the Ensembl database (https://www.ensembl.org/, accessed on 30 May 2022). The genome assembly v1.0 in Ensembl is the same used for data analysis by Anastasiadi and Piferrer [39], however, gene annotation has since been updated according to the Ensembl Gene Annotation pipelines. To facilitate comparative epigenomic analysis with humans, we converted the list of genes with differentially methylated regions (DMRs) to the Ensembl genebuild released version from April 2020. To do this, the genomic coordinates (chromosome, start, end position) of DMRs and surrounding 5000 bp regions were intersected with the genebuild Dicentrarchus_labrax.seabass_V1.0.101.gtf. Chromosome names were as in the primary assembly. A total of 1181 unique genes with DMRs were identified in early domesticates.

Anatomically Modern Humans (AMH)
A detailed map of the evolutionary dynamics of DNA methylation in human groups was recently published [34]. DMRs specific to the AMH-lineage as compared to other hominin lineages, i.e., Denisovan and Neanderthal, were identified using a conservative approach to minimize false positives and variability due to factors such as sex or age, as well as DNA methylation data from chimpanzee samples. AMH-lineage DMRs are a set of 873 DMRs that overlap with the gene body or the promoter up to 5000 bp upstream of 588 genes (Supplementary Dataset S2 of Gokhman et al. [34]; Figure S1). The list of genes with DMRs was supplied by [34] with UCSC (University of California Santa Cruz) identifiers (IDs) and we used the https://biotools.fr/human/ucsc_id_converter (accessed on 30 May 2022) tool to convert them to Ensembl IDs to facilitate comparative epigenomics with the European sea bass.

Neurodevelopmental Cognitive Disorders
WS has a clear genetic origin with the hemideletion of 28 genes at 7q11.23. Some of these genes, e.g., BAZ1B, are involved in epigenetic regulation, such as chromatin remodeling, providing a link to the impact on epigenomic patterns in WS [13]. Differential DNA methylation between patients with WS and healthy individuals as controls has been reported in at least two cases in the literature. DMRs identified using the Infinium ® HumanMethylation450 BeadChip array (Illumina, San Diego, CA, USA) in the blood of 20 WS patients vs. 15 healthy controls found DMRs intersecting 551 unique genes [54]. Differentially methylated cytosines (DMCs) were detected more recently in the blood of a larger sample of 90 WS patients vs. 34 healthy controls using the same array and these intersected with 143 unique genes [55]. The two gene lists were combined for further analysis as genes differentially methylated (DM) in WS, with a total of 624 different genes.
SZ is a complex psychiatric disorder and epigenome-wide association studies (EWAS) have been carried out to explore the role of DNA methylation in SZ pathophysiology, with discordant results. Recently, a meta-analysis of five EWAS datasets was published, including samples taken from different parts of the brain (frontal cortex, cerebellum, hippocampus, and prefrontal cortex), between 3 and 47 samples per study and using either the Infinium HumanMethylation450 BeadChip or HumanMethylation27 BeadChip (Illumina, San Diego, CA, USA) [56]. A total of 513 genes were commonly DM in combinations of 4-5 EWAS and these were used here for further analysis as the DM genes in SZ.
ASD refers to a group of complex neurodevelopmental disorders with heterogeneous symptoms and underlying etiology. ASD heritability is complex and genetic variants involved are diverse with their number ranging between 1000 and 3000 genes reflecting ASD heterogeneity [57]. Other molecular aspects to better understand ASD include epigenetic variants and several studies were published in the last years. This allowed us to apply more stringent criteria for inclusion in this study, mainly a minimum number of 15 samples and identification of DMRs, which are considered more robust than DMCs only. Genes from four studies published in the last 4 years, thus, included: (a) 31 genes with DM that were at the same time differentially expressed and common in three independent studies based on blood samples [58], (b) 181 core genes with DMRs detected using all three approaches in blood cells [57], (c) 145 unique differentially expressed genes with DMRs in blood cells from three ASD subphenotypes (severe, intermediate, mild) and a group of combined cases [59], and (d) 58 genes with DMRs detected in postmortem brain samples [60]. The four datasets combined led to a list of 411 unique ASD genes.

Comparative Analyses
The BioMart data mining tool from Ensembl was used to identify orthologues of human genes from the genome assembly GRCh38.p13 in the European sea bass genome. Duplicate entries were eliminated for further analysis. Thus, we identified unique orthologues as follows: 589 for AMH, 506 for WS, 532 for SZ, and 367 for ASD ( Figure S1). The BioMart tool was used to identify paralogues of the human genes involved in neurodevelopmental cognitive disorders in the human genome (GRCh38.p13), in turn used to identify orthologues in the European sea bass genome. Duplicate entries from the combined list of original orthologues and orthologues of human paralogues were eliminated and the number of homologs finally available for comparative analyses was as follows: 3460 for WS, 4000 for SZ and 2994 for ASD.
Pairwise comparisons were performed with the fish early domesticates (FED) as a reference and one human group as its pair. Thus, four pairwise comparisons occurred every time: (1) FED vs. AMH, (2) FED vs. WS, (3) FED vs. SZ, and (4) FED vs. ASD. Overlaps between gene lists were identified and visualized using the InteractiVenn tool [61]. The significance of overlap was tested using Fisher's exact test for testing the independence of two variables represented by a contingency table. As the genomic background for gene overlap testing, the total number of 23,382 genes in the European sea bass genome (Ensembl genebuild released April 2020) was set. Furthermore, we performed Monte Carlo permutations to test whether overlaps were higher than expected by chance. Random samples of genes were drawn without replacement from the 23,882 total gene list according to the specific gene list each time, e.g., to test the overlap of orthologues FED vs. AMH, 1181 genes for FED vs. 589 genes for AMH were randomly drawn in each iteration. The process was repeated 10,000 times and each time the length of the intersection or overlap between the two genes lists was counted. The standard score of permutation was calculated as observed-mean(permuted)/sd(permuted) and the p-value as: times permuted overlap is higher than observed overlap divided by the number of permutations (10,000). Fisher's tests and permutations were performed using R (v. 4.0.0) [62] and Rstudio (v. 1.4.1717) [63].
The Enrichr tool was used for enrichment analyses and knowledge discovery of gene sets [64][65][66]. Enrichment analyses were performed for the initial lists of genes (FED, AMH, WS, SZ and ASD). Enriched pathways from the databases BioPlanet, Wikipathway, Mammalian Phenotype and GO-terms Biological Process were kept for further comparisons, which included overlap testing as previously carried out, with background on the total number of terms found in each library on Erichr. Reduction and visualization of GO-terms were aided by REViGO [67]. IDs of pathways were entered in InteractiVenn to detect overlaps and Fisher's exact tests were run to detect the statistical significance of the overlap. Enrichment analyses were also performed for the genes that overlapped in a pairwise manner between FED genes and homologs (combined lists of orthologues and orthologues of paralogues).

Differentially Methylated Genes during Early Domestication in European Sea Bass and in Humans Are Shared
The early stages of domestication are expected to be associated with DNA methylation changes. To compare DNA methylation changes associated with the early stages of domestication between European sea bass and human, two gene lists were retrieved. In FED 1181 genes with DMRs were detected as compared to wild fish. For humans, based on limited availability and accessibility to early AMH domesticate samples, DNA methylation patterns of present-day AMHs compared to other hominins and primates were considered the most relevant proxy. A total of 589 genes with DMRs were detected as orthologues of AMH. We detected an overlap of 45 genes between FED and AMH and this was significant (Fisher's test, odds ratio = 1.577, p = 0.004; Figure 1a). Furthermore, we found 1.7 times more genes in common between the two gene lists than expected by chance alone (z-score = 13.62, p = 3 × 10 −4 ; Figure 1b). Among the genes with DMRs in both groups (Table 1), we detected several genes that were repeatedly found to be involved in domestication in several species. For example, ADAM metallopeptidases with thrombospondin type 1 motifs, ephrin (eph) receptors, members of the integrin family (alpha or beta), or fibroblast growth factor receptors have been detected in other domesticates (see Dataset 1 from [39] Anastasiadi and Piferrer Among the genes with DMRs in both groups (Table 1), we detected several genes that were repeatedly found to be involved in domestication in several species. For example, ADAM metallopeptidases with thrombospondin type 1 motifs, ephrin (eph) receptors, members of the integrin family (alpha or beta), or fibroblast growth factor receptors have been detected in other domesticates (see Dataset [68][69][70][71][72] for each species). One of these genes is nuclear factor I X (NFIX in humans and nfxib in fish) which was found to be in the top 10 genes with DMRs in AMHs showing a strong correlation between methylation and expression [34]. Several lines of evidence suggest that hypermethylation of NFIX associates with its downregulation in the AMH lineage [34]. In FEDs, nfixb was hypermethylated in the testis (+29.98%) but hypomethylated in the muscle tissue (−35.88%). In other tissues, other nuclear factor 1 isoforms contained DMRs: in muscle tissue, nuclear factor 1 a-type contained 2 DMRs with opposite methylation patterns (+20.69% and −42.30%), and in brain tissue, nuclear factor 1 a-type contained 2 hypomethylated DMRs (−27.36% and −34.02%) and nuclear factor 1 b-type isoform X2 contained a hypomethylated DMR (−30.14%).

Early Domestication in European Sea Bass and Neurodevelopmental Cognitive Disorders Affect Paralogue Genes
Genes exhibiting DNA methylation changes in patients with neurodevelopmental cognitive disorders with traits parallel to the domestication syndrome such as SZ, WS and ASD were obtained from the literature. A total of 532 genes with DM were orthologues to SZ patients, 506 genes with DM in WS patients and 367 genes with DM in ASD patients.

Early Domestication in European Sea Bass and Neurodevelopmental Cognitive Disorders Affect Paralogue Genes
Genes exhibiting DNA methylation changes in patients with neurodevelopmental cognitive disorders with traits parallel to the domestication syndrome such as SZ, WS and ASD were obtained from the literature. A total of 532 genes with DM were orthologues to SZ patients, 506 genes with DM in WS patients and 367 genes with DM in ASD patients. These gene lists of orthologues were compared to the genes of FED to evaluate whether DNA methylation in common genes was affected by these conditions. The pairwise overlaps were not significant in all cases, with 28 genes overlapping in SZ (odds ratio = 1.04, p = 0.439; Figure S2a), 31 overlapping in WS (odds ratio = 1.233, p = 0.88; Figure S2b) and 23 genes overlapping in ASD (odds ratio = 1.262, p = 0.169; Figure S2c). Permutation testing for the pairwise comparisons showed that the number of overlaps was within the range expected by chance in the case of SZ (z-score = −0.59, p = 0.216; Figure S2d) and ASD (z-score = 2.61, p = 0.062; Figure S2f), and only marginally significant in the case of WS (z-score = 3.75, p = 0.049; Figure S2e).
In an attempt to overcome the constraints of the conservative approach applied here for orthologues and since key candidate genes of domestication were present in all pairwise comparisons, e.g., protocadherins, ADAM metallopeptidases, collagens and glutamate receptors, we then focused on comparisons of functional properties. Orthologue genes were submitted for enrichment analyses and pairwise comparisons were performed at the pathway level following the reasoning that similar processes may be affected by different genes. We considered four libraries targeted by Enrichr as the most informative in our case: Bioplanet, WikiPathways, GO-terms Biological Process, and MGI Mammalian Phenotype. Terms in all four libraries were examined for enrichment according to the gene lists we provided (FED, SZ, WS, and ASD) and pairwise comparisons of terms were performed as follows: (1) FED vs. SZ, (2) FED vs. WS, and (3) FED vs. ASD ( Figure S3). In 42% of the comparisons, there was no overlap of terms, while in three cases there were between 1 and 4 terms overlapping. The overlaps of terms were significant only in cases of SZ for WikiPathways (odds ratio = 4.477, p = 0.003; Figure S3d) and GO Biological Process (odds ratio = 2.442, p = 0.002; Figure S3g). WikiPathways included endochondral ossification with skeletal dysplasia (WP4808) and endochondral ossification (WP474) like in the enrichment of orthologue genes overlapping in AMH, but also neural crest differentiation (WP2064). GO Biological Process enriched included development of the renal system (GO:0072001), kidney (GO:0001822), or ureteric bud (GO:0001657), as well as regulation of immune cells such as T-helper 17 and alpha-beta T (GO:2000317, GO:0046639, or GO:2000320). Taken together these results indicate that further comparative analyses could reveal more additional similarities.
To investigate the role of gene families, we compared gene lists containing not only the orthologues but also the paralogues of genes. The FED gene list was maintained in the original format and served as the control in the pairwise comparisons completed as above. For the other three gene lists (SZ, WS, and ASD), paralogues in the human genome were obtained by Biomart, merged with the original genes and then orthologues in the European sea bass genome were identified, resulting in lists containing unique homologues (orthologues and paralogues). The gene lists contained 4000 homologues for SZ, 3460 homologues for WS, and 2994 homologues for ASD. Overlap between all pairwise comparisons was significant with 241 genes common in SZ (odds ratio = 1.258, p = 0.001; Figure 3a, Dataset S1), 236 in WS (odds ratio = 1.470, p = 4.422 × 10 −7 ; Figure 3b, Dataset S2), and 178 overlapping in ASD (odds ratio = 1.222, p = 0.011; Figure 3c, Dataset S3). Since these gene lists contain~8 times more genes than previously, the significance of the overlaps could be attributed to larger numbers. To test whether the number of overlaps could be expected by chance due to large number of genes, we performed Monte Carlo permutations using random sampling of genes from the whole genome as previously. We found that overlaps between gene lists were higher than expected by chance in all cases, including SZ (z-score = 49.03, p = 0; Figure 3d  To evaluate the functional properties of the core overlaps between genes in FED and lists of homologous genes in cognitive disorders, we performed enrichment analysis using Enrichr as previously. Pathways affected in all pairwise comparisons included neural crest differentiation (WP2064), ectoderm differentiation (WP2858), hair follicle development: organogenesis-part 2 of 3 (WP2839), arrhytmogenic right ventricular cardiomyopathy (WP2118; Figure 4a-c, full lists in Tables S4-S6). Pathways affected in at least two pairwise comparisons included endochondral ossification with skeletal dysplasia (WP4808) like in the core overlap of FED with orthologues of AMH, or also focal adhesion (WP306) and BMP signaling in eyelid development (WP3927) among others (Figure 4a-c).   Tables S7-S9). In SZ and WS, the extracellular matrix organization was the most significantly enriched GO-term. In ASD, the most significantly enriched GO-term was renal system development and among the enriched GO-terms, we detected glutamatergic synaptic transmission (Figure 5c), a process involving glutamate receptors that have been recognized as affected by domestication across species [39,53].  Tables S7-S9). In SZ and WS, the extracellular matrix organization was the most significantly enriched GO-term. In ASD, the most significantly enriched GO-term was renal system development and among the enriched GO-terms, we detected glutamatergic synaptic transmission (Figure 5c), a process involving glutamate receptors that have been recognized as affected by domestication across species [39,53].

Discussion
We have shown that a sizeable portion of epigenetic changes in early European sea bass domesticates occur in similar genes when compared to AMHs, and in similar gene families as in human-specific neurodevelopmental cognitive disorders. Thus, parallel epigenetic changes seem to manifest in independent processes across vertebrates involving domestication, hypothesized to be self-domestication in the case of humans. Our finding that similar genes or gene families exhibited epigenetic changes between human groups and European sea bass provides evidence for domestication as a process affecting similar functional biological properties in vertebrates. Further, it indicates that fish species can be suitable models for research on epigenomics in the context of HSD, as well as human cognitive disorders.
For the purposes of this study, we implemented in silico analyses to compare the lists of genes that exhibited epigenetic changes, measured as differences in DNA methylation. We followed a very conservative approach and included layers of statistical testing, however, some inevitable limitations associated with the nature of the study are present. Genes with epigenetic changes have been pulled from different studies which have used distinct methodologies to interrogate methylation status (e.g., arrays or sequencing) and distinct algorithms to analyze them. In the case of AMH, the data were deduced from comparisons with reconstructed methylomes using a robust methodology. However, this dataset lacks actual methylation data for early AMHs, as comparisons were performed

Discussion
We have shown that a sizeable portion of epigenetic changes in early European sea bass domesticates occur in similar genes when compared to AMHs, and in similar gene families as in human-specific neurodevelopmental cognitive disorders. Thus, parallel epigenetic changes seem to manifest in independent processes across vertebrates involving domestication, hypothesized to be self-domestication in the case of humans. Our finding that similar genes or gene families exhibited epigenetic changes between human groups and European sea bass provides evidence for domestication as a process affecting similar functional biological properties in vertebrates. Further, it indicates that fish species can be suitable models for research on epigenomics in the context of HSD, as well as human cognitive disorders.
For the purposes of this study, we implemented in silico analyses to compare the lists of genes that exhibited epigenetic changes, measured as differences in DNA methylation. We followed a very conservative approach and included layers of statistical testing, however, some inevitable limitations associated with the nature of the study are present. Genes with epigenetic changes have been pulled from different studies which have used distinct methodologies to interrogate methylation status (e.g., arrays or sequencing) and distinct algorithms to analyze them. In the case of AMH, the data were deduced from comparisons with reconstructed methylomes using a robust methodology. However, this dataset lacks actual methylation data for early AMHs, as comparisons were performed with Neanderthal and Denisovans, which are considered as non-self-domesticated hominin species [34]. With regards to neurodevelopmental diseases, due to their often complex etiology, there may be differences in genes detected as DM by different authors who may have used different sampling strategies. Thus, we chose to include only studies which fulfilled stringent criteria. For example, studies involving a very small number of samples, such as comparisons of a pair of twins, were excluded. However, we cannot rule out the possibility that the exact gene lists with epigenetic changes may vary slightly when following consistent and unified guidelines for their detection. Furthermore, to detect homologs, the Biomart tool from the Ensembl database was used, which is one of the most transparent approaches to perform the task since versions of the genome and annotations can be traced. For the enrichment analyses, genes have to be well characterized and included in the query databases to be informative of the affected pathways. Relying on these bioinformatics resources carries the inherent risks of minor modifications in future updated versions. Nevertheless, the results of this study can be interpreted while taking these limitations into account, since they are based on conservative inclusion criteria and statistical testing and can be used as a step for further research on comparative epigenomics between phylogenetically distant vertebrates.
As noted in the Introduction, the HSD hypothesis, as well as the involvement of neural crest cells in HSD, even though attractive, remained mostly supported theoretically until recently. Genomic approaches comparing genes under positive selection between domesticated mammals and AMHs are starting to be used as supportive evidence for the HSD hypothesis [12,73]. Recently, the hypothesis was empirically validated and the role of BAZ1B, a gene within the hemideleted region in WS with an established role in neural crest induction and migration, was demonstrated [13]. The implication of this gene in morphological and behavioral phenotypes typical of the domestication syndrome via neural crest cell development was further shown using zebrafish as a model [74]. Our comparative results between AMHs and early European sea bass domesticates provide additional support for the role of specific genes in key processes with an impact on (self-)domestication features and suggest a role for epigenetic regulation of their expression. That said, one limitation of our approach is that comparative (epi)genomic approaches should be taken with more caution in comparison to conclusive mechanistic studies. However, appropriate model systems for experimental studies are impossible to obtain for the HSD, as viable cells from early AMH are unavailable, and extremely difficult to obtain for other domesticated vertebrates. Thus, in silico comparative studies, such as this one, may provide insightful information on the genes undergoing molecular changes. An additional problem arises from the fact that the number of genes may increase with the number of phenotypic traits to be considered, which is large in the case of domestication and the domestication syndrome, and which, up to a degree, is species-specific [3]. Greater phylogenetic distance, like between humans and fish, may explain more dissimilarities between traits and, thus, more genes altered overall but with lower biological significance individually to explain domestication. These limitations need to be taken into account when interpreting comparative analyses.
With all these limitations in mind, specific genes involved in key processes have been underlined by our analyses. NFIX is associated with craniofacial skeletal disease phenotypes and related to speech capabilities, and it has already been highlighted for its role in the development of the AMH face and larynx [13]. Another gene common between early domesticated European sea bass and AMH was GLI family zinc finger 3 (GLI3) which is a known transcriptional repressor involved in tissue development, including limb development, and immune system development [75]. GLI3 has a role during embryogenesis, controlling thalamic development [76], as well as calvarial suture development [77], while in ≈98% of Altaic Neanderthals and Denisovans it contains a mutation that is mildly disruptive [78]. The RUNX family transcription factor 3, RUNX3, is involved in the developing spinal cord and also has a role in the language and social regions of brain [79,80]. SMOC1, as well as SMOC2, play a role in endochondral bone formation and are regulated by another member of the RUNX family transcription factor [81]. RUNX2 encodes a master transcription factor during vertebrate development involved in the globularization of the human skull/brain. RUNX2 is also involved in the development of thalamus, which is functionally connected to many genes that are important for brain and language development, and that have experienced changes in our recent evolutionary history [74]. NCOR2 has already been identified as under selection in dogs [69], and is part of the cranial neural crest gene expression program [82]. The above-mentioned genes participate in the enriched mammalian phenotypes detected which match the domestication syndrome traits, like abnormal cranium morphology, hypopigmentation or decreased body strength, but also in human-distinctive features potentially associated with our self-domestication. Similarly, GLI3 and SMOC1 participate in the enriched GO-terms processes related to limb development, including limb morphogenesis and embryonic digit morphogenesis. The GO-term most significantly enriched according to its p-value ranking was the negative regulation of alpha-beta T cell differentiation. This is likely due to the involvement of the above-mentioned genes, i.e., RUNX3, GLI3 and SMOC1, in the immune system as well. These results together reinforce the role of epigenetics in the regulation of similar genes associated with the domestication syndrome during the early stages of domestication in the absence of deliberate selection, as is the case in both humans and fish. These results also provide support for the view that domestication constitutes an example of "developmental bias", i.e., when perturbed by an altered environment, complex organisms pursue a limited number of developmental pathways [4].
As also highlighted in the Introduction, neurodevelopmental cognitive disorders in humans have been previously suggested as models for testing the HSD hypothesis [9,19,20]. WS has already been used to gather molecular evidence for the shaping of the human face and behavior underlying HSD [13]. Our initial analyses in search of common genes and pathways epigenetically altered in European sea bass early domesticates and cognitive disorders were unsuccessful. However, even though orthologue genes seemed to be absent, it was evident that similar gene families were affected, thus, justifying our subsequent approach in the search of paralogues. The lack of common genes could be due to the phylogenetic distance between species and to the nature of conditions tested, i.e., disease phenotypes vs. fish under farming conditions. In SZ and WS, genes of key families were affected including, ADAM metallopeptidases, bone morphogenetic proteins, ephrins, fibroblast growth factors, homeoboxes, laminins, and members of the TBC1 domain family. ADAM metallopeptidases and laminins constitute the core members of the most significantly enriched GO-term of overlapping genes in both comparisons: extracellular structure organization. The role of DM genes of the extracellular matrix has already been highlighted in relation to early domestication in fish, and especially for DM changes established already early during development [39]. At the same time, the brain extracellular matrix is known to have multiple roles in brain development and function, and abnormal alteration of this matrix is increasingly acknowledged as a key etiological factor involved in neurological and psychiatric disorders (see Dityatev et al. [83] for review). In ASD, genes were slightly different and included bone morphogenetic proteins, glutamate receptors, laminins, protocadherins, and semaphorins. Neural crest migration depends on the interaction of receptors, e.g., ephrins and receptors for bone morphogenetic proteins, with extracellular matrix molecules, e.g., laminins and semaphorins [84]. The term "neural crest differentiation" was enriched in the overlapping groups of genes and consistently found in all three neurodevelopmental disorders, together with ectoderm differentiation, hair follicle development: organogenesis-part 2 of 3, and arrhytmogenic right ventricular cardiomyopathy. Members of this term were fgfr2, pax3, axin2, hdac10, cdh2, hes1, tfap2a, tfap2b, and tcf7l1. Disorders of the processes related to the neural crest are often regarded as underlying SZ, WS, and ASD [84]. FGF has an essential embryonic function during vertebrate development and Fgf signalling has been shown to serve as a target for selection during the domestication [85]. In ASD, paralogues of two key genes found in the AMH comparison were also identified as epigenetically altered, i.e., runx3 and gli3. This reinforces the idea that parallel processes are involved in HSD phenotype emergence, either evolutionary or pathologically, supporting the view that cognitive diseases can result from changes in genes involved in the human evolution [17,18]. Together these results show that epigenetic changes occur in similar gene families in independent models of early (self-)domestication and that several of these genes have already an established role in the neural crest and other processes recognized as affected by (self-)domestication.
Fish as animal models have long been used in basic science. Small teleost fish, like zebrafish or medaka, has been recently considered as models to study human neurological disorders including ASD [86], peripheral neuropathy [87], and behavioral neuroscience [88] since they possess several key advantages [89]. First, they consist of a phylogenetically diverse group with species that have evolved phenotypes naturally mimicking human diseases, called "evolutionary mutant models" [90][91][92]. Cross-species comparisons allow for the identification of the best models to study a specific physiological pathway [43]. Furthermore, in model species like zebrafish, genetic mutants for specific genes can be easily generated. Second, since they are vertebrates, their brain's basic structure and function exhibit similarities to humans showing conserved neuronal circuitry [93]. Third, teleost genomes show homology with 70% of genes associated with human diseases [94,95]. Fourth, model fish species larvae are transparent, offering the opportunity for direct observation of the central nervous system during the development [96]. Thus, the use of fish models to study neurodevelopmental cognitive disorders exhibiting (self-)domestication-related features already has a sound basis in previous research. Indeed, zebrafish has been used as a model for the three disorders studied here, SCZ [97], WS [98], and ASD [99]. Further, fish species have been suggested as models to investigate evolutionary questions [100], and their potential as models for domestication has been recently recognized [101]. Our findings that homolog genes were differentially methylated in both human disorders and early European sea bass domesticates provide further evidence for the use of fish species as models to study the epigenomic regulation implicated in HSD-related phenotypes, which has proven to be a key source of the human uniqueness [6].
For research related to the HSD hypothesis, fish not only possess the above-mentioned advantages but also show a key similarity distinct from most farm animals: fish domestication and HSD took place in absence of deliberate selection. Our result that DNA methylation changes in European sea bass early domesticates and human groups manifested in overlapping genes supports the implication of epigenetic mechanisms in domestication as a process of adaptation to a human-made environment, including the environment resulting from our self-domestication. A recent study using zebrafish investigated the role of neural crest in the morphological and behavioral domesticated phenotypes in [13] HSD. They found that a loss of function of the key gene in WS and for the neural crest, baz1b, identified as important previously in humans as well [13], resulted in mild neural crest deficiencies during development and behavioral changes related to stress and sociality in adulthood [74]. Furthermore, comparative genomics using domesticated mammals have already been used to shed light to the HSD hypothesis [12]. Together these results show that fish species can be implemented in comparative (epi)genomics approaches and functional studies to further shed light on the validity of the HSD hypothesis.

Conclusions
We have demonstrated the occurrence of parallel epigenetic changes during independent domestication events in phylogenetically distant vertebrates. These events were driven by living in human-made environments, including the creation of the very humanspecific niche through self-domestication, rather than by intentional selection. Epigenetic changes could be the first level of response to a new environment and could later be genomically integrated. An important part of these parallel epigenetic changes arises in genes associated with the neural crest, further supporting the involvement of mild deficits during neural crest development in the emergence of the domestication syndrome. Other common epigenetic changes manifest in genes with neurological or morphological functions that have been associated with the domestication phenotype, including HSD. These findings contribute to our understanding of the initial molecular changes happening during early (self-)domestication and pave the way for future studies using fish as models to investigate epigenetic changes as drivers of HSD, but also as etiological factors of human-specific cognitive diseases.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13060987/s1, Figure S1: Conceptual design of the study; Figure S2: Overlap of orthologue genes differentially methylated in fish early domesticates (FED) and cognitive disorders; Figure S3: Overlap of pathways enriched associated with orthologue genes differentially methylated in fish early domesticates (FED) and cognitive disorders; Table S1: Enrichment of Gene Ontology (GO) Biological Process (BP) terms associated with genes shared between early fish domesticates and anatomically modern humans; Table S2: Enrichment of Mammalian Phenotype (2014) terms associated with genes shared between early fish domesticates and anatomically modern humans; Table S3: Enrichment of Wikipathways associated with genes shared between early fish domesticates and anatomically modern humans; Table S4: Enrichment of Wikipathways associated with genes shared between early fish domesticates and human groups with schizophrenia; Table S5: Enrichment of Wikipathways associated with genes shared between early fish domesticates and human groups with Williams syndrome; Table S6: Enrichment of Wikipathways associated with genes shared between early fish domesticates and human groups with autism spectrum disorders; Table S7: Enrichment of Gene Ontology (GO) Biological Process (BP) terms associated with genes shared between early fish domesticates and human groups with schizophrenia; Table S8: Enrichment of Gene Ontology (GO) Biological Process (BP) terms associated with genes shared between early fish domesticates and human groups with Williams syndrome; Table S9: Enrichment of Gene Ontology (GO) Biological Process (BP) terms associated with genes shared between early fish domesticates and human groups with autism spectrum disorders; Dataset S1: homologue genes with epigenetic changes in human groups with schrizophrenia; Dataset S2: homologue genes with epigenetic changes in human groups with Williams syndrome; Dataset S3: homologue genes with epigenetic changes in human groups with autism spectrum disorders.