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

DNA methylation is one of the key regulators of eukaryotic genomes, and can both inhibit [ 1 , 2 ] and enhance gene expression [ 3 , 4 ], depending on where the DNA methylation occurs. This DNA methylation can be environmentally driven [ 5 ], but can also be modified and regulated via DNA variation [ 6 , 7 , 8 , 9 , 10 ]. We have previously addressed this using a wild x domestic chicken ( Gallus gallus domesticus ) 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 [ 4 ]. However, the corresponding knowledge of 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 and sex balancing on the Z chromosome in the chicken, this is particularly relevant. Sex chromosome dosage compensation [ 11 , 12 ] is a process that prevents the expression imbalance originating from the differential number of sex chromosomes present in males or females when homo- and hetero-gametic sexes exist, relative to the autosomes. 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 in some ways distinct from sex balancing, which is to balance the expression of genes on the sex chromosomes between males and females. However, the two processes may share at least a partial mechanism in common, in particular in instances where dosage compensation acts to decrease gene expression in the homogametic sex on the sex chromosomes. Although this may harmonise gene expression between the sexes on the sex chromosomes, it also results in half the expression of the sex chromosomes relative to the autosomes. In the case of XX/ XY systems, these have chromosome wide regulation (see [ 13 ]), though the mechanism between species is not always the same. In the case of Drosophila , up-regulation of the X in males (XY) occurs, resulting in similar expression between males and females for the sex chromosome(s), and also sex chromosomes and autosomes having similar expression, regardless of sex [ 14 ]. In the case of placental mammals (also XX/ XY), the X chromosome in the homogametic sex (females) is reduced via X chromosome inactivation [ 15 , 16 , 17 ]. This results in both sexes having a uniform ratio of expression for the sex chromosomes, and also that autosomes have on average twice the level of expression relative to the sex chromosomes for both sexes (Pessia 2014). This is achieved via epigenetic mechanisms, notably DNA methylation and histone modification [ 18 ]. Dosage compensation is less well-described in ZW systems, with the female typically being heterogametic, in contrast to the mammalian XY systems [ 15 , 16 ]. ZZ/ ZW systems lack full global Z chromosome compensation, meaning that the homogametic sex has increased expression of the sex chromosome, relative to the heterogametic sex [ 19 ]. Very incomplete and location-specific dosage compensation does appear to occur [ 20 , 21 , 22 ]. In birds, gene expression on the Z in males (ZZ) is not double that of females, but instead genes on the Z are on average around 30% upregulated in males at the transcript level [ 20 , 21 , 23 ]. Similarly, males (ZZ), have a similar sex chromosome expression relative to the autosomes, whilst in females this is reduced. In other ZZ/ ZW systems an effect more similar to XX/ XY systems is observed. In the silk moth, B . mori , the homogametic sex has reduced expression, bringing them in-line with the heterogametic sex chromosome expression, thereby reducing both sexes sex chromosome expression relative to the autosomes [ 12 ]. 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 [ 24 ], whilst more recently an additional region on 73.16–73.17Mb was also identified [ 25 ]. 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. [ 26 ], two pooled samples from Whole Genome Bisulfite sequenced chicken were used in [ 27 ], and one male and one female White Throated Swallow was used in [ 25 ]). This makes it harder to detect smaller regions, and in particular the scope of inter-individual 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 [ 28 ]. Large-scale analysis of within species variation could give a better resolution of hypermethylated regions as well as detect differences between individuals in sex-specific 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 [ 6 ], 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 [ 4 ]. 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 sex-balancing related dosage compensation, 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 DNA methylation is used to regulate sex-differences in gene expression on the Z chromosome. In particular, we are interested in how DNA methylation balances the expression of genes on the sex chromosome between males and females (primarily, rather than the balancing of expression between sex chromosomes and autosomes). For this purpose, an intercross between lines that strongly diverge for DNA methylation levels is an excellent model, as it creates a population with a wide degree of variation in both DNA methylation and gene expression. In this case, the origin of the methylation differences (wild vs domestic) is not as important for sex balancing aspects as the wide variation in DNA methylation present, enabling correlations between widely varying DNA methylation and gene expression levels on the Z chromosome to be assessed in both sexes. We hypothesise that DNA methylation plays a distinct role in sex balancing in the MHM region, beyond its role in dosage compensation. Thus, whilst the MHM region is known for the almost complete suppression of nine lncRNAs in males, it also has the capacity to quantitatively modify Z chromosome gene expression variation. If this is distinct to dosage compensation, then potentially sex balancing can decrease the expression of the homogametic sex. In this case, whilst genes on the Z may be balanced between sexes, they are not compensated between the autosomes and the Z chromosomes in either sex. Further, we hypothesise that the regulation of sex chromosomes by DNA methylation in birds is at least partly regional, as opposed to the global DNA methylation-mediated mechanism described in mammals. One potential benefit of some regional mechanisms is that it allows for the resolution of intra-locus sexual conflict on the sex chromosome, by means of trans-regulation. The role that DNA methylation plays in sex balancing on the Z chromosome can be investigated by the relative abundance of sex interaction effects in cis and trans methylation QTL that are either present on the Z chromosome (cis), or present on the autosomes but affect DNA methylation on the Z chromosome (trans). Similarly, sex interaction effects in cis and trans expression QTL affecting the Z chromosome can also be investigated. In these cases, sex balancing can be far more easily achieved via sex interaction effects, although these can also interact with dosage compensation effects. We would thus predict that sex interaction effects affecting the Z chromosome should be highly prevalent, though if these effects interact with Z chromosome dosage compensation, then local methylation and expression QTL may possess fewer sex interactions as they are reduced by the compensation effects. Similarly, we would also predict that where sex balancing occurs, this is due to decreased expression of the homogametic sex, rather than increased expression of the heterogametic sex, if sex balancing and dosage compensation are distinct to one another.

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 S1 Table ). The original MHM region was hypothesised to be approximately 460kb in length [ 24 ]. 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 Fig 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). Male methylation is shown in blue, female methylation is shown in red. In total, 19 MHM regions (hereafter referred to as blocks) were identified (see Table 1 and Figs 2 and S1 ). Of these continuous blocks, 17 had genes in the local vicinity. In this instance, we defined local as being within 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. To test if dosage compensation acts locally on a gene-by-gene basis or uniformly throughout each block, the methylation levels within these MHM blocks were correlated with the neighbouring genes (see Table 1 ), i.e. individual methylation windows present 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 Figs 2 and S1 . Interestingly, neighbouring genes frequently displayed different correlations with methylation (i.e. neighbouring genes could have very different correlations with local DNA methylation), indicating that these regions seem to be associated with expression on a gene-by-gene basis. 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. 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. 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 Fig 2 ) that are almost entirely silent in males, this region (see Figs 2 and S1 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 Fig 2 , and MHM#12,13,14,15,16,19 in S1 Fig ). Therefore, increased methylation in males is associated with a reduction in the differences in male-biased gene expression, but does 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. One other MHM region has previously been putatively identified at 73.16–73.173Mb on the Z chromosome by Sun et al. [ 25 ]. 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-difference, but the median level of methylation over all individuals was lower than in the original MHM region). Nevertheless, the region shows very significant DNA methylation levels differences between the sexes (see S2 Table ), with significantly more male DNA methylation. All of the neighbouring genes to these MHM regions were also assessed for potential GO enrichment, with no GO enrichment found for those genes in the immediate vicinity. 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 Tables 2 and S3 ). 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. 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. To test if females up-regulate gene expression in the MHM regions, or if males down-regulate gene expression, genes that were present in the combined MHM regions were divided into those that had a significant sex difference and those that did not (see Table 1 ). Fourteen genes had a significant sex difference, and 22 did not (though as noted above, many still had less than the average male: female expression ratio). For the genes that displayed significant sex differences, males had an average expression of 5952 (S.D. 8804), whereas for balanced genes (no sex difference) the average male expression was 343 (S.D. 323). A similar pattern was seen for females, with unbalanced genes having an average expression of 4869 (S.D. 7227), and balanced genes having an average expression of 308 (S.D. 276). Of the two genes that displayed female biased expression (i.e. increased expression in females), the average male expression was 335 (S.D. 318), and the average female expression was 1643 (S.D. 1977). For both males and females, balanced genes had a significantly lower expression (males, t-test statistic -2.4, p = 0.03, females t-statistic = -2.4, p = 0.035). This potentially indicates that balancing occurs by decreasing male gene expression, though female gene expression is also low for these genes. Note that one issue here is that the unbalanced genes in the MHM region are still being balanced to some extent (with lower gene expression differences on average than the Z chromosome on average). To assess the role that methylation plays in sex balanced and sex unbalanced genes in the MHM regions, we tested for a difference in the number of correlations for each gene with the local methylation windows present in their respective MHM regions. In the case of unbalanced genes within the MHM regions, an average of 4.64 correlations were found per gene (S.D. 2.4). All of these correlations were positive (i.e. methylation was positively correlated with gene expression). For the unbalanced female-biased genes within the MHM region, 9 correlations were found per gene (S.D.0), with all of these having a negative correlation with methylation. In the case of the balanced genes, 0.95 correlations per gene were found (S.D. 1.9). Of these, 0.18 correlations per gene had a negative correlation and 0.77 correlations per gene had a positive correlation. Therefore sex-balanced genes (no significant sex difference in gene expression levels) show significantly fewer correlations with local MHM-region methylation (t-statistic = -4.8, p = 6.7x10 -5 ), and whilst unbalanced genes all show positive correlations between methylation and gene expression, balanced genes display a mixture of positive and negative correlations. For female biased unbalanced genes, this group display the most correlations per gene, with all of these being negative. These two genes are the lncRNAs that were previously identified in the initial MHM region. These results appear to show that quantitative changes in methylation are only correlated with gene expression when no (or less) sex balancing occurs, with increasing local methylation appearing to increase, rather than decrease, gene expression. When genes are sex-balanced, local methylation appears to play less of a role in regulating gene expression. 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 S4 Table . 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 S5 Table . 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 [ 4 ], 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 S4 Table ) and 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. 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. 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 [email protected] contained genes enriched for immunoglobulin-fold/domain, the hotspot (ID#5 in Table 3 ) on [email protected] had genes enriched for the activity of glutathione and metabolism of cytochrome P450, and the hotspot (ID#6 in Table 3 ) on [email protected] contained genes enriched for activity with rhodopsin, see S6 Table . The hotspots and their distribution across the genome are illustrated in Fig 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). (A) The 12 autosomal hotspots affecting Z DNA methylation, and (B) the single Z chromosome hotspot affecting Z and autosomal methylation. To test the role that cis and trans regulation plays in sex-balancing gene expression on the Z chromosome, we selected all the cis-acting methQTL that were present on the Z chromosome, as well as any trans acting loci that affected methylation on the Z chromosome. These QTL were then divided into those with a significant sex-interaction, and those without one. In this case, a sex interaction indicates that the allelic effect of the QTL is different depending on the sex of the recipient, i.e. an interaction between QTL and sex is identified. Note that this is distinct from the fixed factor effect of sex (i.e. when one sex always has a uniformly higher or lower expression relative to the other). Thirty-eight cis methQTL were present on the Z chromosome, and none of these had a significant sex-interaction. In stark contrast, 540 trans methQTL affected Z chromosome methylation, and 528 of these displayed a significant sex interaction, with only 22 not displaying such an interaction (chi-squared test, p<0.0001). This ratio is similar even if only significant and not suggestive loci are considered. In this case, 45 trans methQTL showed a sex interaction effect, and only 3 did not. To test how sex balanced and sex unbalanced genes are differently regulated over the whole Z chromosome, all of the 51 eQTL that were either present on the Z chromosome or affected gene expression on the Z chromosome (but were trans in effect) were assessed. The trans effects that were located on the autosomes but affected gene expression on the sex chromosome (n = 43) all had a significant sex interaction. On average, gene expression for these genes in males was 3276 (S.D. 6521), and 2923 (S.D. 6082) in females. In stark contrast, all the cis-acting eQTL that were present on the Z chromosome had no significant sex interaction effect. Similarly, all the trans-acting eQTL that were based on the Z chromosome, but affected autosomal gene expression, also had no significant sex interaction. The average gene expression values of the genes with no significant sex interaction was 381 (S.D. 549) for males, and 243 (S.D. 246) for females. 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 [ 37 ]. 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 [ 36 ], to assess the orientation of the observed correlation. In this case, this means testing whether the causality relationship indicates that methylation is driving gene expression, or the reverse (or if the two are unrelated). Four of the overlaps had a LEO.NB.OCA-score > 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 , chrZ : 44706094–44707218 ) 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. 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.

Sex balancing is not only vital to regulate basic gene expression between the sexes on the sex chromosomes, but also as a mechanism to resolve intra-locus sexual conflict. Our results indicate that the MHM region contains several examples of quantitative sex-balancing effects occurring in the genes present there. This implies that this region is not purely limited to a means of removing expression of key lncRNAs in males. We also find that sex-balancing is not limited to just the MHM region, but also occurs throughout the Z chromosome via trans-effect autosomally-based QTL. These results agree with the study by Catalan et al [ 77 ] that found methylation was spread throughout the Z chromosome, and that questioned the role of the MHM in dosage compensation. Whilst our results don’t address dosage compensation, they do indicate that sex balancing also operates on a gene-by-gene basis, at least to some extent, and that the autosomes may have an important role in such regulation. In contrast, Z chromosome-based modifier loci appear to be limited in their ability to modify Z chromosome gene expression independently between the sexes.

The study was approved by the Local Ethical Committee (Linköping) of the Swedish National Board of Laboratory Animals. Dnr 16110–2020. 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 [ 29 ]. 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 [ 30 ]. 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 [ 4 ]. 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 [ 30 ] 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 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 [ 31 ]. 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 [ 32 ]. Details of the genetic marker locations can be found in Johnsson et al [ 33 ]. Interval mapping was performed using the “qtl2” R-package [ 34 ]. 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. To prevent QTL that were strongly dependent on outliers, we firstly removed all QTL that had fewer than 10 individuals in any genotype class, and also removed QTL that were driven by large effect single outliers via visual inspection. Gene expression QTL (eQTL) analysis was performed using R/qtl, using RMA preprocessed [ 35 ] 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 [ 30 ], 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. The MHM region was identified using the transcript deposited in the NCBI GenBank by [ 24 ], 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). 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 [ 36 ]. 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. NEO assesses causality by using the SNP markers from each QTL as anchors that are then used for orientating the edges of the network and integrating gene expression and DNA methylation, through the use of Structural Equation Model comparisons (SEM). The extent of edge orientation is calculated with a Local SEM-based Edge Orienting score (LEO), that is compared to the next best model (NB). 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. Direction of causality was tested as follows: eQTL marker -> gene expression -> DNA methylation gene expression <- DNA methylation gene expression -> DNA methylation gene expression <- DNA methylation <- methQTL marker gene expression DNA methylation gene expression -> DNA methylation gene expression <- DNA methylation

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. (PDF) 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. (XLSX) 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. (XLSX) The three FHM blocks (continuous regions) are highlighted in orange and indicated with their block ID in a separate column. 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. (XLSX) The phenotype of each methylation QTL (methylation window), nearest marker to the methylation QTL, LOD score, confidence interval, and nearest marker to each confidence are given, as well as whether the QTL is cis or trans in effect, are all given. (XLSX) Closest marker, LOD score, confidence interval, presence or absence of sex interaction, and nearest marker to the confidence interval are all presented. (XLSX) Category, GO term, p-value (absolute and also FDR controlled), genes involved and fold enrichment are all given. (XLSX)