The regulation of methylation on the Z chromosome and the 2 identification of multiple novel Male Hyper-Methylated regions 3 in the chicken.

22 DNA methylation is a key regulator of eukaryote genomes, and is of particular 23 relevance in the regulation of gene expression on the sex chromosomes, with a key 24 role in dosage compensation in mammalian XY systems. In the case of birds, dosage 25 compensation is largely absent, with it being restricted to two small Male Hyper-26 Methylated (MHM) regions on the Z chromosome. To investigate how variation in 27 DNA methylation is regulated on the Z chromosome we utilised a wild x domestic 28 advanced intercross in the chicken, with both hypothalamic methylomes and 29 transcriptomes assayed in 124 individuals. The relatively large numbers of individuals 30 allowed us to identify additional genomic MHM regions on the Z chromosome that 31 were significantly differentially methylated between the sexes. These regions appear 32 to down-regulate local gene expression in males, but not remove it entirely (unlike 33 the lncRNAs identified in the initial MHM regions). In addition, trans effect hotspots 34 were also identified that were based on the autosomes but affected the Z, and also 35 that were based on the Z chromosome but that affected autosomal DNA methylation


Introduction
DNA methylation is one of the key regulators of eukaryotic genomes, and can both inhibit (Gaston andFried 1995, Mann, Chatterjee et al. 2013) and enhance gene expression (Yin, Morgunova et al. 2017, Höglund, Henriksen et al. 2020), depending on where the DNA methylation occurs.This DNA methylation can be environmentally driven (Nalvarte, Rüegg et al. 2018), but can also be modified and regulated via DNA variation (Kasowski, Kyriazopoulou-Panagiotopoulou et al. 2013, Kilpinen, Waszak et al. 2013, McVicker, van de Geijn et al. 2013, Bélteky, Agnvall et al. 2018, Guerrero-Bosagna 2019).We have previously addressed this using a wild x domestic chicken model to study the regulation of variation in autosomal DNA methylation, and how it can quantitatively regulate gene expression using a QTL mapping based approach.This enabled us to identify how domestication in the chicken led to a small number of large-effect trans hotspots, where these loci regulated variation in DNA methylation throughout the genome.Moreover, we found methylation to not only be the driver but also the response to gene expression variation (Höglund, Henriksen et al. 2020).However, the corresponding regulation of DNA methylation variation on the Z chromosome is still lacking.For example, the extent to which quantitative variation in DNA methylation is controlled between the autosomes and sex chromosomes is an open question, as is the extent to which DNA methylation is regulated on the Z chromosome in general.Given the role that DNA methylation plays in dosage compensation on the Z chromosome in the chicken, this is particularly relevant.Dosage compensation prevents the expression imbalance originating from the number of sex chromosomes present in males or females when homo-and heterogametic sexes exist.Dosage compensation occurs when the dose effect due to one sex having only a single sex chromosome, and therefore half the number of gene copies, is compensated by either decreasing gene expression in the homogamete or increasing expression in the heterogamete.This can be over the whole sex chromosome or over specific regions (Mank 2013).Dosage compensation is less well-described in ZW systems, with the female typically being heterogametic, in contrast to the mammalian XY systems (Itoh, Melamed et al. 2007, Vicoso andBachtrog 2011).In the mammalian XY system, dosage compensation is achieved by X inactivation, achieved via epigenetic mechanisms, notably DNA methylation and histone modification (Fang, Disteche et al. 2019).However, such chromosomal inactivation is largely absent from birds, with instead very incomplete and locationspecific dosage compensation, if any (Ellegren, Hultin-Rosenberg et al. 2007, Mank andEllegren 2009).Despite this, gene expression on the Z in males is not double that of females, but instead genes on the Z are on average around 30% upregulated in males (Ellegren, Hultin-Rosenberg et al. 2007, Melamed and Arnold 2007, Mank and Ellegren 2009).
DNA methylation still plays an important role for sex difference regulation on the avian Z chromosome.In particular, the Male Hyper-Methylated (MHM) region at 27.3Mb was first discovered by Teranishi and colleagues (Teranishi, Shimada et al. 2001), whilst more recently an additional region on 73.16-73.17Mbwas also identified on the Z chromosome (Sun, Maney et al. 2019).With the initial region, it was found that males had greatly increased methylation in an approximately 500kb area, with nine genes that were present there not being expressed in males.In the case of the more recently discovered MHM region at 73.16Mb (designated MHM2), this was smaller and contained three long non coding RNAs (lncRNAs) that were female-biased in expression.In general, these studies are based on small numbers of samples, generally focussing on between species comparisons (for example, one great tit sample was used in Laine et al. 2016, two pooled samples from Whole Genome Bisulfite sequenced chicken were used in (Zhang, Yan et al. 2017), and one male and one female White Throated Swallow was used in (Sun, Maney et al. 2019)).
This makes it harder to detect smaller regions, and in particular the scope of interindividual variation in these MHM regions.This is concerning, particularly considering the degree of DNA methylation variation across individuals in populations and the role of methylation in phenotype formation (Heyn, Moran et al. 2013).Large-scale analysis of within species variation could give a better resolution of hypermethylated regions as well as detect differences between individuals in sexspecific methylation and gene regulation.Various questions still remain regarding the MHM regions, and the genes they contain.The sizes of the MHM regions and the effects of the decreased gene expression is particularly noteworthy -are these genes involved in fundamental sex differences?Similarly, are the genes within these MHM regions regulated in a region-by-region basis or on a gene-by-gene basis?Gene expression regulation via methylation is not restricted to solely promoter regions (Kasowski, Kyriazopoulou-Panagiotopoulou et al. 2013), but can affect gene expression (both positively and negatively) due to effects at enhancer sites, Transcriptional Elements (TEs) and the like.For example, our previous study based on autosomal methylation variation in the chicken found that there was a bias to being positively correlated, whilst correlations between methylation and gene expression could be found within a megabase upstream and downstream of the gene itself (Höglund, Henriksen et al. 2020).Given this, how far away from these MHM regions are genes being affected?Is this still affecting dosage compensation if it upregulates genes?To investigate how DNA methylation variation is regulated on the Z chromosome, as well as the potential role of methylation in dosage compensation and sex differences, we conducted a DNA methylation quantitative trait locus (methylation QTL) analysis using an advanced intercross between domestic chickens and wild Red Junglefowl.We assayed the hypothalamic transcriptome and methylome on the Z chromosome for 124 individuals, having previously assayed the autosomes for these individuals.It was therefore possible to map both cis and trans related loci that modulate variation in DNA methylation on the Z chromosome, as well as to assess how methylation is used to regulate sex-differences in gene expression on the Z chromosome.

Methods
The study population was composed of 124 chickens (55 females, 69 males) from which the hypothalamus tissue was dissected out at day 212.The individuals used were from an 8th generation advanced intercross, founded using a Red Junglefowl (wild) male and three White Leghorn (domestic) females.A detailed description of the intercross generation, housing conditions, etc can be found in (Johnsson, Gustafson et al. 2012).

RNA and DNA methylation isolation
RNA was isolated from the hypothalamus tissue which was homogenised using Ambion TRI Reagent (Life Technologies) following the manufacturer's protocol.cDNA synthesis and microarray-based gene expression were performed using a Nimblegen 135k array, as described previously (Johnsson, Williams et al. 2016).DNA was isolated from the remainder of the TRI reagent homogenate by mixing 125µl ice-cold 99% ethanol with 250µl TRI reagent homogenate.Samples were vortexed, incubated on ice for 5min and centrifuged at 12'000 RPM for 10 min in room temperature.The pellet was saved and isolation continued using the DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer's protocol.DNA methylation was assessed by Methylated DNA immunoprecipitation (MeDIP) protocol.Further details of the MeDIP protocol can be found in (Höglund, Henriksen et al. 2020).

Phenotypes: methylation and gene expression
DNA methylation phenotypes were generated by dividing the chicken genome into 1000bp windows, yielding a total of 1050176 methylation windows, of which 82426 were located on chromosome Z.The MeDIP-seq reads were mapped to each methylation window and normalised by dividing with the total read count for each individual respectively.Sequencing was performed on an IonProton machine (Thermo Fisher Scientific) using the Torrent Suite software (version 4.4.1) by the National Genomics Infrastructure in Uppsala, Sweden.The sequence depth was on average 3.4X ± 0.97 (standard deviation), the read length was on average 136 ± 15 bp, the raw reads was on average 23.8 million ± 5.2 and the quality score was on average 22 ± 1.The Gene expression dataset has been published previously (Johnsson, Williams et al. 2016) and was based on the NimbleGen 12 x 135K Custom Gene expresson array, mapping to 22628 unique genes composed of Ensembl, RefSeq genes and Expressed Sequence Tags.

Quantitative Trait Loci (QTL) analysis
Quantitative Trait Loci (QTL) analysis was performed to identify genomic regions associated with the variation found within DNA methylation levels for the 1 million methylation windows.A genetic marker map was generated using 652 SNP markers, of which 542 were fully informative between the original parental animals used to generate the intercross.Average marker distance was ~16 cM, as per recommendations (Darvasi and Soller 1994).Of these, 36 markers were present on the Z chromosome with a 15cM average marker distance.Note that as the intercross is a linkage-based cross and not a GWAS of an outbred population (which relies on linkage disequilibrium and has built up historical recombinations over hundreds of generations) far fewer markers are required to cover the genome, as it is only required to identify the recombinations that have accrued during the intercrossing (Lynch and Walsh 1998).Details of the genetic marker locations can be found in Johnsson et al (Johnsson, Rubin et al. 2014).Interval mapping was performed using the "qtl2" R-package (Broman, Gatti et al. 2019).This package was used as it is able to correctly analyse sex chromosomes in an advanced intercross.A local (cis) scan was performed, restricted to methylation windows present on the Z chromosome, with the local region considered to be within 50cM up-and down-stream of each methylation window.A trans scan was also performed.In the case of the trans scan, a full genome scan was performed for trans effect methylation QTL that were located on either the autosomes or Z chromosome that affected methylation on the Z chromosome.In addition, a scan was also performed for trans methylation QTL located on the Z chromosome that were associated with methylation present on the autosomes.Sex and batch were set as covariates in the test model, with sex also used as an interactive covariate, where significant (if the LOD score of the sex interaction model was >1 LOD higher than the non-sex interaction model).
Significance thresholds were determined via a permutation test with and without sex interactions for both local (cis) and trans methylation QTL.Local (putatively cis) regions were defined as 50cM up and downstream to the closest genetic marker, whilst anything outside this region was defined as trans.For the trans permutations, 20000 random methylation phenotypes were permuted 1000 times each, both for sex and non-sex interaction, and for cis permutations 17000 random phenotypes permuted 1000 times each.From the permutations the top 5 % LOD-scores for each phenotype were saved and from these the top 5% were chosen as significance threshold and the top 20% as the suggestive threshold, respectively.This yielded significance cis LOD-score of 5.73 (sex interaction) and 4.29, (no sex interaction), with suggestive thresholds of 4.87 and 3.58.For the trans thresholds, significance was at LOD-score of 7.70 and 7.68 (sex and non-sex interaction, respectively), whilst the suggestive threshold was 5.92 and 5.93.
Gene expression QTL (eQTL) analysis was performed using R/qtl, using RMA preprocessed (Irizarry, Bolstad et al. 2003) expression levels as quantitative phenotypes with sex and batch as additive covariates.The same criteria for cis-eQTL was applied as for the autosomes (see (Johnsson, Williams et al. 2016), with local eQTL defined as those within +/-50cM of the gene, with trans referring to any other location.Significance thresholds for cis and trans eQTL were 4.0 and 6.0, respectively.

Male Hyper-Methylated (MHM) region
The MHM region was identified using the transcript deposited in the NCBI GenBank by (Teranishi, Shimada et al. 2001), accession AB046698 (2332 bp), with this being the probe sequence used to identify the region initially.This sequence maps to two genomic locations: chrZ:27375241-27391116 (99.1% match) and chrZ:27329191-27333743 (98.9% match), hereafter referred to as MHMa and MHMb respectively.
The MHMa and MHMb regions were corroborated in our dataset and the parameters for methylation levels obtained were used to identify other MHM-like regions.These parameters were: median methylation status per window of > 8.52, a sex difference equal to a Wilcoxon rank sum test/Mann-Whitney test p-value < 1.75e-10 and comprising of five or more adjacent methylation windows (i.e.these values are those identified for the original MHM region in our dataset).

QTL overlaps
With both mQTL and eQTL mapped it was possible to assess whether any correlations could be found between DNA methylation levels and gene expression which are both associated to a locus.By overlapping the confidence intervals of the mQTL and eQTL, and regressing the gene expression with methylation, genomic regions that putatively control either the methylation or gene expression (or both) were observed.The correlation was tested with all individuals and sex as a factor, and with the sexes separate, yielding 3 models.Any genes that significantly correlated with a methylation window were finally tested for causality using the Network Edge Orientation (NEO) package in R (Aten, Fuller et al. 2008).In this way, it is possible to ascribe hypothetical orientation of the regulatory relationship, whether DNA methylation regulates gene expression or vice versa.Significance using the NEO package is based on the LEO.NB score, which quantifies the support of the best fitting causal model versus the second best fitting model.As both the eQTL and methylation QTL originated from the same genotype (imputed marker position) and thus are treated as a single-marker orientation with a LEO.NB.OCA-score > 1.0 considered significant, and a score of > 0.8 as suggestive.

Dosage Compensation and the Male Hyper Methylated (MHM) Region
To assess the degree of male-biased hyper methylated regions, we first analysed the previously known hyper methylated regions -MHMa and MHMb.These two regions, situated at 27.375Mb and 27.329Mb respectively, had a 3.3 and 3.6-fold increase in methylation in males, respectively, with these ratios being highly significant (max Wilcoxon pvalue < 4e-10 and 1.7e-10, respectively, for each region), see Supplementary Table 1).The original MHM region was hypothesised to be approximately 460kb in length (Teranishi et al. 2001).When we assessed the methylation around these two regions, we find elevated methylation from 27.142Mb-27.40Mb(259-kb long), more accurately demarking this region, see Figure 1.To identify further male biased methylation windows, we performed a chromosome-wide scan calculating the degree of sex bias.Based on the pre-existing MHM region, we then selected all those regions with both a strongly significant sex bias (p<1.75e-10, as compared to the average sex bias in methylation on the Z chromosome being p=0.14 and a 1.69 fold methylation difference between males and females) and with at least five adjacent methylation windows (see methods section).
In total, 19 MHM regions (hereafter referred to as blocks) were identified (see Table 1, Figure 2 and Supplementary Figure 1).Of these continuous blocks, 17 had genes in the local vicinity.In this instance, we defined local as being with 100kb of the MHM block, as in our previous study we found strong correlations between gene expression and DNA methylation even up to 100kb away from the gene itself.within each block were correlated with the expression of adjacent genes, controlling for multiple testing.Of the 17 blocks with adjacent genes, 14 had a significant correlation between at least one methylation window and local gene expression, see Figure 2 and Supplementary Figure 1.Interestingly, neighbouring genes frequently displayed differential correlation with methylation, indicating that these regions seem to be associated with expression on a gene-by-gene basis.In total, 51 unique genes (38 present in our dataset) were adjacent to these MHM-like blocks, with 224 significant correlations with methylation levels (methylation windows) of which 134 correlations were negative and 90 positive (tvalue from linear model).Furthermore, of the 38 genes present in our dataset, 34 had a significant sex bias expression with 20 being expressed higher in males and 14 higher in females (M:F ratio).The average fold difference between males and females on the Z chromosome was 1.22 while for the autosomes this was 1.02.In the case of the original MHM region, apart from the RNAse genes (EST probes X603141644 and X603862378 for the lncRNA ENSGALG00000051419 in Figure 2) that are almost entirely silent in males, this region (see Figure 2, Supplementary Figure 1, and Table 1) also contains multiple genes that are still male-biased, but below the average degree of male-bias on the Z chromosome.Similarly, these genes tend to be positively correlated with local methylation, where such a correlation exists.This pattern is also replicated in the newly identified MHM regions (see MHM#1 and #2 in Figure 2, and MHM#12,13,14,15,16,19 in Supplementary Figure 1).Therefore, increased methylation in males is associated with a reduction in the differences in male-biased gene expression, but not eliminate it entirely, in both the existing and the new MHM regions.None of the methylation QTL detected on the Z chromosome (either QTL or phenotypes) overlapped with these MHM regions.

(TABLE 1 PRESENTED SEPARATELY)
One other MHM region has previously been putatively identified at 73.  on the Z chromosome by Sun et al. (2019).We also identify this region in our data, though the median methylation threshold fell slightly below the threshold we set, and was therefore excluded initially (i.e.there was a strongly significant sex-327 difference, but the median level of methylation over all individuals was lower than in 328 the original MHM region).Nevertheless, the region shows very significant DNA 329 methylation levels differences between the sexes (see Supplementary Table 2), with 330 significantly more male DNA methylation.All of the neighbouring genes to these 331 MHM regions were also assessed for potential GO enrichment, with no GO enrichment found for those genes in the immediate vicinity.

Female Hyper-Methylated Regions
As well as additional MHM regions, a search for regions with a lower than average male: female methylation ratio was also performed to identify regions that showed a relative decrease in DNA methylation in males or an increase in DNA methylation in females.Using a criterion of a significant increase in female methylation, relative to males, we firstly identified a total of 118 1kb windows that were significantly more methylated in females than males (see Table 2 and Supplementary Table 3).Of these, three regions consisted of five or more consecutive female-biased methylated windows.These regions were located at 30195000-30200000bp, 42633000-42638000bp, and 49073000-49073000bp on the Z chromosome.No genes were found in these regions, however.An overlap between methylation QTL and these regions was also performed, though once again no overlaps occurred.

Methylation QTL present on the Z chromosome
Methylation QTL were assessed by performing local (cis) methylation QTL scans restricted solely to the Z chromosome.In addition, trans scans were also performed, where the QTL was located on the Z chromosome, but the target methylation window was free to be present on either the Z chromosome or the autosomes.In total, we identify 18 significant cis methylation QTL and 53 significant trans methylation QTL that are based on the Z chromosome, with a further 20 suggestive cis methylation QTL and 528 trans methylation QTL.As expected, most of the methylation QTL (n=528) had a significant sex interaction effect.This is expected due to the large differences in Z chromosome methylation between males and females, with males possessing two methylated chromosomes (ZZ) and females only one (ZW).A full list of all methylation QTL can be found in Supplementary Table 4.In addition, 51 expression QTL (eQTL) were identified on the Z chromosome (either as a QTL or the trans-effect phenotype of a QTL), see Supplementary Table 5. ready.

Trans Methylation QTL Hotspots Affecting the Z chromosome
To identify trans-acting hotspots, we identified where multiple methylation QTL were associated with the same marker and had overlapping confidence intervals.Of the 619 methylation QTL on the Z chromosome, these mapped to 141 different SNP loci.Of these loci, 13 were associated with multiple methylation windows/phenotypes (10 or more methylation windows associated with each marker, respectively).These hotspots on average spanned 5.87Mb of physical distance in the genome (found by taking the shared overlapping confidence intervals and finding the minimum overlapping size), see Table 3.Of note, all bar one (n=12) of these trans hotspots were located on the autosomes, but regulated variation in methylation on the Z chromosome.Of these 12, 3 were previously identified as regulating methylation variation on the autosomes in this intercross (Höglund et al. 2020), on chromosomes 3 (at 18Mb, hotspot 4), 6 (at 7.7Mb, hotspot 6) and 7 (at 2.4 Mb, hotspot 9).One hotspot was located on the Z chromosome (at 41.7Mb, hotspot 13, with this hotspot spread over three adjacent SNPs, rs16768340, rs16782623, rs14016786, see Supplementary Table 4) regulated variation in methylation on different windows in the Z chromosome, as well as some methylation windows on the autosomes.Thus, whilst the majority of regulation in methylation variation appears to be located on the autosomes, with these loci then regulating methylation on the Z chromosome, there is also some regulation of methylation variation by the Z chromosome itself, and even a small amount of autosomal regulation from the Z chromosome.The genes present within these hotspots were further checked for potential enrichment via gene ontology analysis, using DAVID 6.8 (https://david.ncifcrf.gov/).In total 3 hotspots showed enrichment using the DAVID 6.8 database: the hotspot (ID#2) at chr1@91.7MB contained genes enriched for immunoglobulin-fold/domain, the hotspot (ID#5 in Table 3) on chr3@86.5Mbhad genes enriched for the activity of glutathione and metabolism of cytochrome P450, and the hotspot (ID#6 in table 3) on chr4@1.3Mbcontained genes enriched for activity with rhodopsin, see Supplementary Table 6.The hotspots and their distribution across the genome are illustrated in Figure 3. Gene enrichment analysis was also performed for the target genomic regions in the vicinity (±10kb) of each methylation window associated with a methylation QTL hotspot.Some enrichment was found for hotspot ID#4 (located on chromosome 3 at 17.98Mb), however, this result was non-significant (Bonferroni p-value > 0.05).

Causality Analysis Between Methylation and Gene Expression on the Z chromosome
In total 360 overlaps were found between eQTL and methylation QTL.These were methylation and expression QTL where either the QTL or methylation phenotype were located on the Z chromosome.The overlapping phenotypes (gene expression and methylation) were tested for association using a linear model.Of these, 15 overlaps were significant after applying an FDR-based multiple testing corrections.
Eleven of the overlaps were significant (p-value < 0.05, FDR corrected) using all individuals, while 3 were significant (p-value < 0.05, FDR corrected) using only females, and 1 was significant (p-value < 0.05, FDR corrected) using only males, see Table 4.These overlaps contained 5 unique probesets belonging to 2 unique genes and 3 ESTs.The gene LINGO1 (ENSGALG00000002708; chr10:3212741-3290778) is an immunoglobulin domain protein (Yang, Jiang et al. 2022).Immunoglobulin activity was also found in the methylation QTL hotspot on chromosome 3.Additionally, the 15 overlaps were tested with NEO, a network edge orientated method which uses the underlying QTL genotype as anchors for the network (Aten et al., 2008) influences the methylation levels in the region of chrZ:45163000-45165000, see Table 4.This gene has been retired on the GalGal6 genome, with no known function.
In addition one further EST (X603865974) was suggestive (LEO.NB.OCA >0.8), with methylation appearing to drive gene expression in this case.However, as the model p-value was significant, this means that other models (gene expression driving methylation) cannot be ruled out.also occurs.Interestingly, three of the hotspots previously identified as regulating DNA methylation in domestication (primarily via reducing DNA methylation in domestic birds) also appear to regulate DNA methylation on the Z chromosome (Höglund 2020).
As well as the regulation of variation in methylation, we also identified additional Male Hyper-Methylated regions present on the Z chromosome.Unlike the initial MHM region found (Teranishi, Shimada et al. 2001), which identified that the lncRNAs present were completely switched off in males, the regions we identify appear to instead decrease male gene expression, though rather than reduce it entirely, it is instead down-regulated to levels more closely found in females (i.e. reduced male gene expression, relative to female gene expression).This is despite the regions having a similar pattern of sex-differentiated methylation as is seen in the original MHM region.Further, the strength of the methylation differences between sexes was greater in the new regions we identified when compared to the region at 74Mb (although we also identify the 74Mb MHM region as well).These genes thus appear to be linked to sex-based differences between males and females.
No methylation QTL overlap these regions, implying that these regions are not responsible for regulating variation in methylation, which would then fit with these regions instead regulating more basal sex-differences rather than betweenpopulation variation.This idea is reinforced when considering the functions of the genes in these regions.
Of the 18 known genes that are present within the MHM regions, their functions can be broadly divided into learning/ behaviour, bone allocation, development, reproduction, growth/ metabolism and methyl transferase activities.
These tie-in well with the known sex-differences that exist in the chicken.Starting with behaviour, strong behavioural differences exist between males and female chickens (Vallortigara, Cailotto et al. 1990, Nätt, Agnvall et al. 2014, Elfwing, Nätt et al. 2015, Bélteky, Agnvall et al. 2018).In particular, females have decreased anxietyrelated behaviour, though this may be test-dependent (Schutz, Kerje et al. 2002, Campler, Jöngren et al. 2009, Johnsson, Williams et al. 2016, Johnsson, Henriksen et al. 2018, Fogelholm, Inkabi et al. 2019).Of the genes present in the MHM, four are related to behaviour or neurogenesis.The gene SLC1A1 has been shown to play a role in obsessive compulsive disorder and sterotype behaviour (Zike, Chohan et al. 2017, Huang, Liu et al. 2021), as well as schizophrenia susceptibility (Horiuchi, Iida et al. 2012, Li, Su et al. 2020).Anxiety behaviour in chickens has previously been shown to be related to schizophrenia, depression and other mood-based disorders in humans, even sharing some of the same susceptibility loci (Johnsson, Williams et al. 2016, Johnsson, Henriksen et al. 2018).Furthermore, the OCD effects arising from SLC1A1 are stronger in males, so sex-differences in the gene effects have already been demonstrated (Wendland, Moya et al. 2009, Veenstra-VanderWeele, Xu et al. 2012).ZDHHC2I is a major palmitoyl acyltransferase, with decreasing expression leading to increased depression-like behaviours (Gorinski, Bijata et al. 2019).Homer1 also has functions relating to learning and memory (Clifton, Cameron et al. 2017), and also causes susceptibility to Alzheimers (Urdánoz-Casado, Sánchez-Ruiz de Gordoa et al. 2021).In the case of the latter, these effects are strongly sexdependent, only occurring in women.Continuing with bone allocation, female chickens have a complex bone allocation, whereby during egg production the hard outer cortical bone is first mobilised into soft, spongy medullary bone in the centre of the femur, before then being transferred to create the egg shell (one of the major limiting factors in egg production) (Bloom, Domm et al. 1958, Mueller, Schraer et al. 1964).Therefore male and female chickens differ markedly in their bone metabolism -males possess almost no medullary bone, whilst female medullary bone deposition is strongly associated with reproductive output (Johnsson, Gustafson et al. 2012, Johnsson, Rubin et al. 2014, Johnsson, Jonsson et al. 2015).Of the genes in the MHM, Homer1 has numerous beneficial effects in osteoblasts including beta-catenin stabilization (Rybchyn, Brennan-Speranza et al. 2021).FST (foillistatin) is also a powerful regulator of bone metabolism (Gajos-Michniewicz, Piastowska et al. 2010).CER1 has also been found to regulate bone mineral density and be associated with fracture risk.Of note, these effects are found to be strongest in post-menopausal women (Koromila, Dailiana et al. 2012, Koromila, Georgoulias et al. 2013).TLE4 is a critical mediator of osteoblasts and runx2-dependent bone development in the mouse (Shin, Theodorou et al. 2021).Finally, NANS affects skeletal development in zebrafish knock-outs (van Karnebeek, Bonafé et al. 2016).Continuing with reproduction-related genes, the gene FST plays a critical role in mouse uterine receptivity and decidualization (Fullerton, Monsivais et al. 2017), whilst the gene JMY mediates spermatogenesis in mice (Liu, Fan et al. 2020) as well as asymmetric division and cytokinesis in mouse oocytes (Sun, Sun et al. 2011).Finally, CLTA4 is involved in the maintenance of chronic inflammation in endometriosis and infertility (Abramiuk, Bębnowska et al.

2021).
The final category of genes present in the MHM regions affected growth and metabolism, whilst two methyltransferase genes were also present.Large differences in growth and bodyweight exist in the chicken, with males often twice the bodyweight of females.The gene FST, as well as affecting reproduction-related phenotypes, also leads to increased muscle weight in mice when over-expressed (Iyer, Chugh et al. 2021).DMGDH affects body growth through insulin-like growth factor (Baker et al 1993), whilst also affecting selenium status in pregnant women (Mao, Vanderlelie et al. 2016).ABHD3 regulates adipogenic differentiation and lipid metabolism (Linke, Overmyer et al. 2020).Finally, BHMT is a methyltransferase, as is

DMGDH.
In the current study, we have restricted the investigation to a single tissue, albeit repeated in 124 individuals.As such, we are confident that these MHM regions and methylation QTL are present in this tissue type.However the ubiquity of these methylated regions in other tissues must be verified.This opens up the possibility that multiple further MHM regions also exist, but are only present in specific tissue types.This could allow fine-scale regulation of sex differences in a tissue-specific manner.We also assess female-specific hyper-methylated regions, but these were found to be very sparse and had very few genes present, suggesting these are of less importance.In summary, by using a large number of replicates that are assessed for all methylated loci on the Z chromosome, we identify both novel MHM regions in an intra-specific/ inter-population framework, as well as the role that domestication plays in the regulation of the Z chromosome and the genes on it.Figure 2. Four of the 19 Novel MHM regions present on the Z chromosome and their effects on gene expression.Panes illustrate regions 1, 2, 9, 12 (selected as being representative of all the regions).Each pane consists of the following: i) The male:female methylation ratio for the 1kb methylation windows that make up the MHM region (each black dot represents the ratio at one methylation window).The red hashed line at the base indicates the average male:female methylation ratio (~1.7).
ii) Male:female gene expression ratio is indicated by the blue dots, one for each gene in the region, with the ratio shown on the left-side y-axis, and the blue hashed line indicating the average male:female gene expression ratio on the Z chromosome (~1.2).
iii) The number of correlations between each gene and the 1kb methylation windows that make up each MHM.The direction of the correlation (positive or negative) is indicated by the bar being above the line (positive, coloured turquoise) or below the line (negative, coloured purple).The number of correlations is indicated on each bar, whilst each gene name is given on the x-axis.

Table Legends
Table 1.Novel Male Hyper-Methylated (MHM) regions identified in the hypothalamus.The 17 MHM regions containing genes are divided into separate regions, with their location, size, number of probesets present initially given.Also included are the average gene expression values for males and females, the p-value of the sex differences in gene expression, the ratio of male:female gene expression, the number of 1kb windows present within the MHM region that correlate with each gene and the direction of that correlation.Table 2. Novel Female Hyper-Methylated regions identified in the hypothalamus.The position, average and median methylation per window, and the average methylation in males and females per window are all given, as well as the significance of the sexdifference and the average male:female fold ratio.Table 3. List of trans methylation QTL hotspots.Table shows the number of methylation QTL present for each hotspot, its chromosome and base-pair position (nearest marker), and the confidence interval of each hotspot.The number of genes present within the intervals as determined by ensembl.org is also given.Table 4. NEO causality of gene regulation of methylation.The probeset and the methylation window being tested, along with their confidence interval is presented.
In addition, the genotype p-value, the sex p-value (also broken down into male and female), as well as the actual causality statistics (leo.nb.oca and cpa and the model p-value) are all shown.Supplementary Table 1.MHMa and MHMb regions.The position, average and median methylation per window, and the average methylation in males and females per window are all given, as well as the significance of the sex-difference and the average male:female fold ratio.Supplementary Table 2.The previously identified MHM region at 73Mb.The position, average and median methylation per window, and the average methylation in males and females per window are all given, as well as the p-value and significance of the sex-difference and the average male:female fold ratio.Note for the significance of the sex difference, these are classified as non-significant, significant (including a multiple testing correction), and significant at the same level as the original MHM region.
, to assess the orientation of the observed correlation.Four of the overlaps had a LEO.NB.OCAscore > 0.3.Both eQTL and mQTL originated from the same genotype (imputed marker position) and thus are treated as a single-marker orientation where a LEO.NB.OCA-score > 1.0 is significant.Hence, our results indicate that the EST X603598164F1 (gene id:ENSGALG00000050497,

433
Using this wild x domestic paradigm to analyse DNA methylation and its regulation 434 on the Z chromosome in the chicken, we firstly identify over 600 methylation QTL 435 that affect methylation on the Z chromosome.Of these, the majority of trans effect 436 loci are located on the autosomes but affecting the Z chromosome.There were also 437 examples of the reverse, with trans methylation QTL deriving from the Z 438 chromosome but affecting DNA methylation on the autosomes.Furthermore, these 439 trans methylation QTL were concentrated into a small number of hotspots located 440 on the autosomes (n=12), with one hotspot also present on the Z chromosome itself, 441 associated with methylation on the Z chromosome and the autosomes, respectively.442 A total of five genes on the Z chromosome were also candidates for causality 443 between gene expression and methylation, with two passing the network-edge-444 based threshold for significance.Of these, one appears to be a retired gene, whilst 445 the other is an EST of no known function, with the former indicating gene expression 446 affects methylation, whilst the latter indicates that methylation was modifying gene 447 expression.448449The nature of the intercross (a wild bird intercrossed with domestics) allows us to 450 identify consistent differences in methylation that exist between wild and domestic 451 chickens and the regions that associate with and potentially regulate them.With 452 regards to the methylation QTL hotspots identified, it is noteworthy that these are 453 almost all based on the autosomes, with only one situated on the Z chromosome 454 itself.Therefore, the regulation of domestication-based phenotypes with loci present 455 on the Z chromosome appears to generally be autosomally regulated, although the 456 reverse (where autosomal gene expression is regulated by the sex chromosomes)

Figure 1 .
Figure LegendsFigure1.MHMa and MHMb regions (used to identify the original MHM region) and sex differences in methylation levels at these and the surrounding area.Male methylation is shown in blue, female methylation is shown in red.Figure2.Four of the 19 Novel MHM regions present on the Z chromosome and their effects on gene expression.Panes illustrate regions 1, 2, 9, 12 (selected as being representative of all the regions).Each pane consists of the following: i) The male:female methylation ratio for the 1kb methylation windows that make up the MHM region (each black dot represents the ratio at one methylation window).The red hashed line at the base indicates the average male:female methylation ratio (~1.7).ii) Male:female gene expression ratio is indicated by the blue dots, one for each gene in the region, with the ratio shown on the left-side y-axis, and the blue hashed line indicating the average male:female gene expression ratio on the Z chromosome (~1.2).iii)The number of correlations between each gene and the 1kb methylation windows that make up each MHM.The direction of the correlation (positive or negative) is indicated by the bar being above the line (positive, coloured turquoise) or below the line (negative, coloured purple).The number of correlations is indicated on each bar, whilst each gene name is given on the x-axis.Figure 3. Circle plot showing the location of trans methylation QTL hotspots that affect DNA methylation variation on the Z chromosome.(A) The 12 autosomal

Figure 3 .
Figure LegendsFigure1.MHMa and MHMb regions (used to identify the original MHM region) and sex differences in methylation levels at these and the surrounding area.Male methylation is shown in blue, female methylation is shown in red.Figure2.Four of the 19 Novel MHM regions present on the Z chromosome and their effects on gene expression.Panes illustrate regions 1, 2, 9, 12 (selected as being representative of all the regions).Each pane consists of the following: i) The male:female methylation ratio for the 1kb methylation windows that make up the MHM region (each black dot represents the ratio at one methylation window).The red hashed line at the base indicates the average male:female methylation ratio (~1.7).ii) Male:female gene expression ratio is indicated by the blue dots, one for each gene in the region, with the ratio shown on the left-side y-axis, and the blue hashed line indicating the average male:female gene expression ratio on the Z chromosome (~1.2).iii)The number of correlations between each gene and the 1kb methylation windows that make up each MHM.The direction of the correlation (positive or negative) is indicated by the bar being above the line (positive, coloured turquoise) or below the line (negative, coloured purple).The number of correlations is indicated on each bar, whilst each gene name is given on the x-axis.Figure 3. Circle plot showing the location of trans methylation QTL hotspots that affect DNA methylation variation on the Z chromosome.(A) The 12 autosomal binding specificities of human transcription factors."Science 356(6337): eaaj2239.Zhang, M., F.-B. Yan, F. Li, K.-R.Jiang, D.-H.Li, R.-L.Han, Z.-J.Li, R.-R.Jiang, X.-J.Liu and X.-T.Kang (2017)."Genome-wide DNA methylation profiles reveal novel candidate genes associated with meat quality at different age stages in hens."Scientific Reports 7(1): 1-15.Zike, I. D., M. O. Chohan, J. M. Kopelman, E. N. Krasnow, D. Flicker, K. M. Nautiyal, M. Bubser, C. Kellendonk, C. K.Jones, G. Stanwood, K. F. Tanaka, H. Moore, S. E. Ahmari and J. Veenstra-VanderWeele (2017)."OCD candidate gene SLC1A1/EAAT3 impacts basal ganglia-mediated activity and stereotypic behavior."Proc Natl Acad Sci U S A 114(22): 5719-5724.