Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

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

Abstract

DNA methylation is a key regulator of eukaryote genomes, and is of particular relevance in the regulation of gene expression on the sex chromosomes, with a key role in dosage compensation in mammalian XY systems. In the case of birds, dosage compensation is largely absent, with it being restricted to two small Male Hyper-Methylated (MHM) regions on the Z chromosome. To investigate how variation in DNA methylation is regulated on the Z chromosome we utilised a wild x domestic advanced intercross in the chicken, with both hypothalamic methylomes and transcriptomes assayed in 124 individuals. The relatively large numbers of individuals allowed us to identify additional genomic MHM regions on the Z chromosome that were significantly differentially methylated between the sexes. These regions appear to down-regulate local gene expression in males, but not remove it entirely (unlike the lncRNAs identified in the initial MHM regions). These MHM regions were further tested and the most balanced genes appear to show decreased expression in males, whilst methylation appeared to be far more correlated with gene expression in the less balanced, as compared to the most balanced genes. In addition, quantitative trait loci (QTL) that regulate variation in methylation on the Z chromosome, and those loci that regulate methylation on the autosomes that derive from the Z chromosome were mapped. Trans-effect hotspots were also identified that were based on the autosomes but affected the Z, and also one that was based on the Z chromosome but that affected both autosomal and sex chromosome DNA methylation regulation. We show that both cis and trans loci that originate from the Z chromosome never exhibit an interaction with sex, whereas trans loci originating from the autosomes but affecting the Z chromosome always display such an interaction. Our results highlight how additional MHM regions are actually present on the Z chromosome, and they appear to have smaller-scale effects on gene expression in males. Quantitative variation in methylation is also regulated both from the autosomes to the Z chromosome, and from the Z chromosome to the autosomes.

Author summary

DNA methylation is a key regulator of eukaryote genomes, and is of particular relevance in the regulation of gene expression on the sex chromosomes, with a key role in dosage compensation in mammalian XY systems. In the case of birds, dosage compensation is largely absent, with it being restricted to two small Male Hyper-Methylated (MHM) regions on the Z chromosome. We utilised a wild x domestic advanced intercross in the chicken, with both hypothalamic methylomes and transcriptomes assayed in 124 individuals, to investigate the role that methylation plays in regulating gene expression on the Z chromosome. Our results highlight how additional MHM regions are actually present on the Z chromosome, and they appear to have smaller-scale effects on gene expression in males. Quantitative variation in methylation is also regulated both from the autosomes to the Z chromosome, and from the Z chromosome to the autosomes. In addition, these MHM regions were further tested and the most balanced genes appear to show decreased expression in males, whilst methylation appeared to be far more correlated with gene expression in the less balanced, as compared to the most balanced genes.

Citation: Höglund A, Henriksen R, Churcher AM, Guerrero-Bosagna CM, Martinez-Barrio A, Johnsson M, et al. (2024) The regulation of methylation on the Z chromosome and the identification of multiple novel Male Hyper-Methylated regions in the chicken. PLoS Genet 20(3): e1010719. https://doi.org/10.1371/journal.pgen.1010719

Editor: Kees Van Oers, Netherlands Institute of Ecology, NETHERLANDS

Received: March 28, 2023; Accepted: January 31, 2024; Published: March 8, 2024

Copyright: © 2024 Höglund et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Gene expression data (generated with Microarray) for the hypothalamus tissue is available in the ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress) under accession number E-MTAB-3154. The genotypes scored for the QTL analysis is available at Figshare [https://doi.org/10.6084/m9.figshare.12803876]. The DNA methylation data (generated with MeDIP) is available at: https://doi.org/10.6084/m9.figshare.12803873 Finally, the code used to analysye the data and the QTL cross file are available here: https://github.com/flerpan01/DNA-methylation-QTL-mapping-chicken-intercross.

Funding: The project was supported by grants from the European Research Council (Advanced grant GENEWELL to PJ and Consolidator grant FERALGEN 772874 to DW), the Swedish Research Council (VR, grant 2018-04421 to DW), Carl Tryggers Stiftelse (to DW), and the Linköping University Neuro-network (to DW). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

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.

Methods

Ethics statement

The study was approved by the Local Ethical Committee (Linköping) of the Swedish National Board of Laboratory Animals. Dnr 16110–2020.

Study population

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 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 [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].

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 [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

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.

Male Hyper-Methylated (MHM) region

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).

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 [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:

  1. eQTL marker -> gene expression -> DNA methylation <- methQTL marker
  2. eQTL marker -> gene expression <- DNA methylation <- methQTL marker
  3. eQTL marker -> gene expression -> DNA methylation
  4. gene expression <- DNA methylation <- methQTL marker
  5. gene expression <- shared QTL marker -> DNA methylation
  6. gene expression -> DNA methylation <- methQTL marker
  7. eQTL marker -> gene expression <- DNA methylation

Results

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 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).

thumbnail
Download:
Fig 1. 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.

https://doi.org/10.1371/journal.pgen.1010719.g001

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.

thumbnail
Download:
Fig 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.

https://doi.org/10.1371/journal.pgen.1010719.g002

thumbnail
Download:
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.

https://doi.org/10.1371/journal.pgen.1010719.t001

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.

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 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.

thumbnail
Download:
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 sex-difference and the average male:female fold ratio.

https://doi.org/10.1371/journal.pgen.1010719.t002

Sex differences in gene expression in the MHM Regions

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).

Methylation correlations in balanced and unbalanced genes in MHM regions

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 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 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.

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 [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.

thumbnail
Download:
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.

https://doi.org/10.1371/journal.pgen.1010719.t003

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.5Mb had genes enriched for the activity of glutathione and metabolism of cytochrome P450, and the hotspot (ID#6 in Table 3) on chr4@1.3Mb 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).

thumbnail
Download:
Fig 3. Circle plot showing the location of trans methylation QTL hotspots that affect DNA methylation variation on the Z chromosome.

(A) The 12 autosomal hotspots affecting Z DNA methylation, and (B) the single Z chromosome hotspot affecting Z and autosomal methylation.

https://doi.org/10.1371/journal.pgen.1010719.g003

Role of cis and trans methQTL in sex balancing on the Z chromosome

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.

Sex differences in eQTL gene expression and the Z chromosome

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.

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 [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.

thumbnail
Download:
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.

https://doi.org/10.1371/journal.pgen.1010719.t004

Discussion

Using this wild x domestic paradigm to analyse DNA methylation and its regulation on the Z chromosome in the chicken, we firstly identify over 600 methylation QTL that affect methylation on the Z chromosome. Of these, the majority of trans effect loci are located on the autosomes but affecting the Z chromosome. There were also examples of the reverse, with trans methylation QTL deriving from the Z chromosome but affecting DNA methylation on the autosomes. Furthermore, these trans methylation QTL were concentrated into a small number of hotspots located on the autosomes (n = 12), with one hotspot also present on the Z chromosome itself, associated with methylation on the Z chromosome and the autosomes, respectively. A total of five genes on the Z chromosome were also candidates for causality between gene expression and methylation, with two passing the network-edge-based threshold for significance. Of these, one appears to be a retired gene, whilst the other is an EST of no known function, with the former indicating gene expression affects methylation, whilst the latter indicates that methylation was modifying gene expression.

We also assessed the role of methylation and gene expression in the regulation of sex-balancing in both the new and old MHM regions. Sex-balanced genes had decreased expression in the MHM region, relative to the unbalanced/ less balanced genes that were present. This suggests that for these genes, male expression is reduced, with a similar effect to mammalian XY systems. The expression of the balanced genes was less than half the expression of the unbalanced genes, implying that these genes are even more tightly regulated in both sexes. The relationship with DNA methylation in sex balanced and sex unbalanced (or less balanced) genes was also assessed, with a strong difference found with balanced genes displaying far fewer correlations with local methylation than unbalanced genes. In addition, unbalanced genes only displayed positive correlations with local methylation. These indicate that strongly sex-balanced genes are not affected by local methylation, but more weakly balanced or unbalanced genes are more likely to be regulated by local methylation (or at least correlated with it).

The presence of sex interactions in both the methylation and expression QTL that were identified can also provide evidence for the regulation of sex balancing on the Z chromosome. With regard to methylation QTL, the local (cis) methQTL were all strongly balanced, whereas the non-Z chromosome based trans methQTL were all unbalanced in that they displayed strong sex interactions with each QTL. This appears to show that locally-regulated methylation is very stable between the sexes–even when there may be differences in total methylation between the sexes at these loci, these are not interacting with allelic effects and regulation. In contrast, trans-effect loci seem distinct from the sex-balancing effects. Almost all of these are located on the autosomes and show modified Z chromosome methylation levels, implying that they are not part of the sex-balancing mechanisms that exist on the Z chromosome, and could even be a way of introducing more variation between the sexes in terms of methylation regulation on the Z chromosome.

Sex differences in expression QTL on the Z chromosome displayed a similar pattern as methylation QTL, with cis-acting loci on the Z displaying sex balance (or more accurately a lack of any sex interaction with the eQTL), whereas trans-acting loci displayed a large degree of sex-interaction when they originated from the autosomes. Also similar to the genes on the MHM region, both male and female expression values for the sex balanced genes were significantly lower than the unbalanced genes. These results imply that when the cis and trans effects originate from the Z chromosome, whatever factor is balancing these genes remains, and that potentially the genes affected are not related to sex-biased expression traits. In contrast, when the QTL are present on the autosomes, whatever factor that balances expression between the sexes is bypassed, with these genes more likely to be related to sex-biased trait expression due to the large differences that can be induced between the sexes by the alleles in question.

Intra-locus sexual conflict arises from differential selection pressures experienced by males and females, when they share a common genetic basis [38]. This can be resolved when a locus has different effect between the sexes, with this being considered as a separate genetic architecture between the sexes [39], with several examples of phenotypic traits having such a distinct genetic architecture between the sexes [40, 41, 42]. Regulatory evolution is a rapid and efficient solution to intra locus sexual conflict [43], where the regulatory elements of a gene evolve to act in a sex-specific manner. Sex-interacting expression QTL and methylation QTL represent exactly this type of differential genetic architecture between the sexes. The results shown here would therefore suggest that it is possible to induce such separate genetic architecture for genes on the Z chromosome, but only when the QTL are located on the autosomes. When regulatory loci are present on the Z chromosome itself, these do not exhibit this differential architecture, and appear instead to always give equal effects between the sexes. Therefore, the autosomal-derived trans effects appear to drive sex differences (and thus reduce sexual conflict), whereas Z chromosome-based cis and trans effects appear to be unable to escape the dosage compensation and/ or other sex balancing effects that exist on the Z chromosome.

The nature of the intercross (a wild bird intercrossed with domestics) allows us to identify consistent differences in methylation that exist between wild and domestic chickens and the regions that associate with and potentially regulate them. With regards to the methylation QTL hotspots identified, it is noteworthy that these are almost all based on the autosomes, with only one situated on the Z chromosome itself. Therefore, the regulation of domestication-based phenotypes with loci present on the Z chromosome appears to generally be autosomally regulated, although the reverse (where autosomal gene expression is regulated by the sex chromosomes) 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 [4].

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 [24], 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 between-population 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 [10, 44, 45, 46]. In particular, females have decreased anxiety-related behaviour, though this may be test-dependent [30,47,48,49,50]. 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 [51,52], as well as schizophrenia susceptibility [53,54]. 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 [30,47]. Furthermore, the OCD effects arising from SLC1A1 are stronger in males, so sex-differences in the gene effects have already been demonstrated [55,56]. ZDHHC2I is a major palmitoyl acyltransferase, with decreasing expression leading to increased depression-like behaviours [57]. Homer1 also has functions relating to learning and memory [58], and also causes susceptibility to Alzheimers [59]. In the case of the latter, these effects are strongly sex-dependent, 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) [60,61]. 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 [29,33,62]. Of the genes in the MHM, Homer1 has numerous beneficial effects in osteoblasts including beta-catenin stabilization [63]. FST (foillistatin) is also a powerful regulator of bone metabolism [64]. 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 [65,66]. TLE4 is a critical mediator of osteoblasts and runx2-dependent bone development in the mouse [67]. Finally, NANS affects skeletal development in zebrafish knock-outs [68]. Continuing with reproduction-related genes, the gene FST plays a critical role in mouse uterine receptivity and decidualization [69], whilst the gene JMY mediates spermatogenesis in mice [70] as well as asymmetric division and cytokinesis in mouse oocytes [71]. Finally, CLTA4 is involved in the maintenance of chronic inflammation in endometriosis and infertility [72].

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 [73]. DMGDH affects body growth through insulin-like growth factor [74], whilst also affecting selenium status in pregnant women [75]. ABHD3 regulates adipogenic differentiation and lipid metabolism [76]. 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.

The use of eQTL and methQTL studies can help to understand sex differences and sex-balancing as related to the Z chromosome, as well as sexual conflict resolution. Our results indicate that Z chromosome sex-balancing regulation appears to decrease male gene expression for the most tightly regulated genes, whilst the role of local methylation seems to be stronger for more weakly balanced genes or unbalanced genes. Given that methylation QTL and expression QTL studies deal explicitly with within-population quantitative variation, they represent an excellent resource with which to understand the more fine-scale sex-balancing that is occurring, as opposed to the very severe dosage compensation and sex-balancing effects seen in the classical MHM regions.

Conclusions

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.

Supporting information

S1 Fig. All of the Novel MHM regions present on the Z chromosome and their effects on gene expression.

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.

https://doi.org/10.1371/journal.pgen.1010719.s001

(PDF)

S1 Table. 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.

https://doi.org/10.1371/journal.pgen.1010719.s002

(XLSX)

S2 Table. 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.

https://doi.org/10.1371/journal.pgen.1010719.s003

(XLSX)

S3 Table. All female hyper-methylated regions.

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.

https://doi.org/10.1371/journal.pgen.1010719.s004

(XLSX)

S4 Table. List of local (cis) and trans methylation QTL present on the Z chromosome.

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.

https://doi.org/10.1371/journal.pgen.1010719.s005

(XLSX)

S5 Table. Expression QTL (eQTL) present on the Z chromosome.

Closest marker, LOD score, confidence interval, presence or absence of sex interaction, and nearest marker to the confidence interval are all presented.

https://doi.org/10.1371/journal.pgen.1010719.s006

(XLSX)

S6 Table. GO enrichments from methylation QTL hotspots.

Category, GO term, p-value (absolute and also FDR controlled), genes involved and fold enrichment are all given.

https://doi.org/10.1371/journal.pgen.1010719.s007

(XLSX)

Acknowledgments

The research was carried out within the framework of the Swedish Centre of Excellence in Animal Welfare Science, and the Linköping University Neuro-network. SNP genotyping was performed by the NGI Sweden, using the Uppsala Sequencing Center (SNP&SEQ Technology Platform). Computations were enabled by resources provided by the Swedish National Infrastructure for

Computing (SNIC) at the Uppsala SNIC Centre, partially funded by the Swedish Research Council through grant agreement no. 2018–05973. Bioinformatic support was provided by NBIS (National Bioinformatics Infrastructure Sweden).

References

  1. 1. Gaston K, Fried M. CpG methylation has differential effects on the binding of YY1 and ETS proteins to the bi-directional promoter of the Surf-1 and Surf-2 genes. Nucleic acids research. 1995;23(6):901–9. pmid:7731802
  2. 2. Mann IK, Chatterjee R, Zhao J, He X, Weirauch MT, Hughes TR, et al. CG methylated microarrays identify a novel methylated sequence bound by the CEBPB| ATF4 heterodimer that is active in vivo. Genome research. 2013;23(6):988–97. pmid:23590861
  3. 3. Yin Y, Morgunova E, Jolma A, Kaasinen E, Sahu B, Khund-Sayeed S, et al. Impact of cytosine methylation on DNA binding specificities of human transcription factors. Science. 2017;356(6337):eaaj2239. pmid:28473536
  4. 4. Höglund A, Henriksen R, Fogelholm J, Churcher AM, Guerrero-Bosagna CG, Martinez Barrio A, et al. The methylation landscape and its role in domestication and gene regulation in the chicken. Nature Ecology and Evolution. 2020;4(12):1713–24. pmid:32958860
  5. 5. Nalvarte I, Rüegg J, Guerrero-Bosagna C. Intrinsic and Extrinsic Factors That Influence Epigenetics. Epigenetics and Assisted Reproduction: CRC Press; 2018. p. 99–114.
    • 6. Kasowski M, Kyriazopoulou-Panagiotopoulou S, Grubert F, Zaugg JB, Kundaje A, Liu Y, et al. Extensive variation in chromatin states across humans. Science. 2013;342(6159):750–2. pmid:24136358
    • 7. Kilpinen H, Waszak SM, Gschwind AR, Raghav SK, Witwicki RM, Orioli A, et al. Coordinated effects of sequence variation on DNA binding, chromatin structure, and transcription. Science. 2013;342(6159):744–7. pmid:24136355
    • 8. McVicker G, van de Geijn B, Degner JF, Cain CE, Banovich NE, Raj A, et al. Identification of genetic variants that affect histone modifications in human cells. Science. 2013;342(6159):747–9. pmid:24136359
    • 9. Guerrero-Bosagna C, editor From epigenotype to new genotypes: Relevance of epigenetic mechanisms in the emergence of genomic evolutionary novelty. Seminars in cell & developmental biology; 2019: Elsevier.
      • 10. Bélteky J, Agnvall B, Bektic L, Höglund A, Jensen P, Guerrero-Bosagna C. Epigenetics and early domestication: differences in hypothalamic DNA methylation between red junglefowl divergently selected for high or low fear of humans. Genetics Selection Evolution. 2018;50(1):13. pmid:29609558
      • 11. Ohno S. Sex chromosomes and sex-linked genes: Springer Science & Business Media; 2013.
        • 12. Gu L, Walters JR, Knipple DC. Conserved Patterns of Sex Chromosome Dosage Compensation in the Lepidoptera (WZ/ZZ): Insights from a Moth Neo-Z Chromosome. Genome Biol Evol. 2017;9(3):802–16. pmid:28338816
        • 13. Ferrari F, Alekseyenko AA, Park PJ, Kuroda MI. Transcriptional control of a whole chromosome: emerging models for dosage compensation. Nature structural & molecular biology. 2014;21(2):118–25. pmid:24500429
        • 14. Gelbart ME, Kuroda MI. Drosophila dosage compensation: a complex voyage to the X chromosome. 2009.
        • 15. Itoh Y, Melamed E, Yang X, Kampf K, Wang S, Yehya N, et al. Dosage compensation is less effective in birds than in mammals. Journal of biology. 2007;6(1):1–15.
        • 16. Vicoso B, Bachtrog D. Lack of global dosage compensation in Schistosoma mansoni, a female-heterogametic parasite. Genome Biology and Evolution. 2011;3:230–5. pmid:21317157
        • 17. Schulz EG, Heard E. Role and control of X chromosome dosage in mammalian development. Current opinion in genetics & development. 2013;23(2):109–15. pmid:23465885
        • 18. Fang H, Disteche CM, Berletch JB. X inactivation and escape: epigenetic and structural features. Frontiers in cell and developmental biology. 2019;7:219. pmid:31632970
        • 19. Mank JE. The W, X, Y and Z of sex-chromosome dosage compensation. Trends in Genetics. 2009;25(5):226–33. pmid:19359064
        • 20. Ellegren H, Hultin-Rosenberg L, Brunström B, Dencker L, Kultima K, Scholz B. Faced with inequality: chicken do not have a general dosage compensation of sex-linked genes. BMC biology. 2007;5(1):1–12.
        • 21. Mank J, Ellegren H. All dosage compensation is local: gene-by-gene regulation of sex-biased expression on the chicken Z chromosome. Heredity. 2009;102(3):312–20. pmid:18985062
        • 22. Mank JE. Sex chromosome dosage compensation: definitely not for everyone. Trends in genetics. 2013;29(12):677–83. pmid:23953923
        • 23. Melamed E, Arnold AP. Regional differences in dosage compensation on the chicken Z chromosome. Genome biology. 2007;8(9):1–10. pmid:17900367
        • 24. Teranishi M, Shimada Y, Hori T, Nakabayashi O, Kikuchi T, Macleod T, et al. Transcripts of the MHM region on the chicken Z chromosome accumulate as non-coding RNA in the nucleus of female cells adjacent to the DMRT1 locus. Chromosome Research. 2001;9(2):147–65. pmid:11321370
        • 25. Sun D, Maney DL, Layman TS, Chatterjee P, Soojin VY. Regional epigenetic differentiation of the Z Chromosome between sexes in a female heterogametic system. Genome research. 2019;29(10):1673–84. pmid:31548356
        • 26. Laine VN, Gossmann TI, Schachtschneider KM, Garroway CJ, Madsen O, Verhoeven KJF, et al. Evolutionary signals of selection on cognition from the great tit genome and methylome. Nature Communications. 2016;7:10474. pmid:26805030
        • 27. Zhang M, Yan F-B, Li F, Jiang K-R, Li D-H, Han R-L, et al. Genome-wide DNA methylation profiles reveal novel candidate genes associated with meat quality at different age stages in hens. Scientific Reports. 2017;7(1):1–15.
        • 28. Heyn H, Moran S, Hernando-Herraez I, Sayols S, Gomez A, Sandoval J, et al. DNA methylation contributes to natural human variation. Genome research. 2013;23(9):1363–72. pmid:23908385
        • 29. Johnsson M, Gustafson I, Rubin C-J, Sahlqvist A-S, Jonsson KB, Kerje S, et al. A Sexual Ornament in Chickens Is Affected by Pleiotropic Alleles at HAO1 and BMP2, Selected during Domestication. PLoS Genetics. 2012;8(8):e1002914. pmid:22956912
        • 30. Johnsson M, Williams MJ, Jensen P, Wright D. Genetical Genomics of Behavior: A novel chicken genomic model for anxiety behavior. Genetics. 2016;202(1):327–40. pmid:26733665
        • 31. Darvasi A, Soller M. Optimum spacing of genetic markers for determining linkage between marker loci and quantitative trait loci. Theoretical and Applied Genetics. 1994;89(2–3):351–7. pmid:24177853
        • 32. Lynch M, Walsh B. Genetics and Analysis of Quantitative traits. Sunderland, MA: Sinauer Associates; 1998.
          • 33. Johnsson M, Rubin CJ, Höglund A, Sahlqvist AS, Jonsson K, Kerje S, et al. The role of pleiotropy and linkage in genes affecting a sexual ornament and bone allocation in the chicken. Molecular ecology. 2014;23(9):2275–86. pmid:24655072
          • 34. Broman KW, Gatti DM, Simecek P, Furlotte NA, Prins P, Sen Ś, et al. R/qtl2: Software for Mapping Quantitative Trait Loci with High-Dimensional Data and Multiparent Populations. Genetics. 2019;211(2):495–502. pmid:30591514
          • 35. Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP. Summaries of Affymetrix GeneChip probe level data. Nucleic Acids Research. 2003;31(4):e15. pmid:12582260
          • 36. Aten JE, Fuller TF, Lusis AJ, Horvath S. Using genetic markers to orient the edges in quantitative trait networks: the NEO software. BMC Systems Biology. 2008;2(1):34. pmid:18412962
          • 37. Yang H, Jiang L, Zhang Y, Liang X, Tang J, He Q, et al. Anti-LINGO-1 antibody treatment alleviates cognitive deficits and promotes maturation of oligodendrocytes in the hippocampus of APP/PS1 mice. The Journal of comparative neurology. 2022;530(10):1606–21. pmid:35014704
          • 38. Bonduriansky R, Chenoweth SF. Intralocus sexual conflict. Trends in ecology & evolution. 2009;24(5):280–8. pmid:19307043
          • 39. Innocenti P, Morrow EH. The sexually antagonistic genes of Drosophila melanogaster. PLoS biology. 2010;8(3):e1000335. pmid:20305719
          • 40. Randall JC, Winkler TW, Kutalik Z, Berndt SI, Jackson AU, Monda KL, et al. Sex-stratified genome-wide association studies including 270,000 individuals show sexual dimorphism in genetic loci for anthropometric traits. PLoS genetics. 2013;9(6):e1003500. pmid:23754948
          • 41. Dapper AL, Wade MJ. The evolution of sperm competition genes: the effect of mating system on levels of genetic variation within and between species. Evolution. 2016;70(2):502–11. pmid:26748568
          • 42. Karp NA, Mason J, Beaudet AL, Benjamini Y, Bower L, Braun RE, et al. Prevalence of sexual dimorphism in mammalian phenotypic traits. Nature communications. 2017;8(1):15475. pmid:28650954
          • 43. Wright AE, Fumagalli M, Cooney CR, Bloch NI, Vieira FG, Buechel SD, et al. Male-biased gene expression resolves sexual conflict through the evolution of sex-specific genetic architecture. Evolution Letters. 2018;2(2):52–61. pmid:30283664
          • 44. Nätt D, Agnvall B, Jensen P. Large sex differences in chicken behavior and brain gene expression coincide with few differences in promoter DNA-methylation. PLoS One. 2014;9(4):e96376. pmid:24782041
          • 45. Elfwing M, Nätt D, Goerlich-Jansson VC, Persson M, Hjelm J, Jensen P. Early stress causes sex-specific, life-long changes in behaviour, levels of gonadal hormones, and gene expression in chickens. PLoS One. 2015;10(5):e0125808. pmid:25978318
          • 46. Vallortigara G, Cailotto M, Zanforlin M. Sex differences in social reinstatement motivation of the domestic chick (Gallus gallus) revealed by runway tests with social and nonsocial reinforcement. Journal of Comparative Psychology. 1990;104(4):361. pmid:2282785
          • 47. Johnsson M, Henriksen R, Fogelholm J, Höglund A, Jensen P, Wright D. Genetics and Genomics of Social Behavior in a Chicken Model. Genetics. 2018;209(1):209–21. pmid:29531010
          • 48. Fogelholm J, Inkabi S, Höglund A, Abbey-Lee R, Johnsson M, Jensen P, et al. Genetical Genomics of Tonic Immobility in the Chicken. Genes. 2019;10(5):341.
          • 49. Campler M, Jöngren M, Jensen P. Fearfulness in red junglefowl and domesticated White Leghorn chickens. Behavioural Processes. 2009;81(1):39–43. pmid:19154782
          • 50. Schutz K, Kerje S, Carlborg O, Jacobsson L, Andersson L, Jensen P. QTL analysis of a red junglefowl x White Leghorn intercross reveals trade-off in resource allocation between behavior and production traits. Behavior Genetics. 2002;32(6):423–33. pmid:12467340
          • 51. Zike ID, Chohan MO, Kopelman JM, Krasnow EN, Flicker D, Nautiyal KM, et al. OCD candidate gene SLC1A1/EAAT3 impacts basal ganglia-mediated activity and stereotypic behavior. Proc Natl Acad Sci U S A. 2017;114(22):5719–24. pmid:28507136
          • 52. Huang X, Liu J, Cong J, Zhang X. Association Between the SLC1A1 Glutamate Transporter Gene and Obsessive-Compulsive Disorder in the Chinese Han Population. Neuropsychiatric disease and treatment. 2021;17:347–54. pmid:33574671
          • 53. Li W, Su X, Chen T, Li Z, Yang Y, Zhang L, et al. Solute Carrier Family 1 (SLC1A1) Contributes to Susceptibility and Psychopathology Symptoms of Schizophrenia in the Han Chinese Population. Front Psychiatry. 2020;11:559210. pmid:33173509
          • 54. Horiuchi Y, Iida S, Koga M, Ishiguro H, Iijima Y, Inada T, et al. Association of SNPs linked to increased expression of SLC1A1 with schizophrenia. American journal of medical genetics Part B, Neuropsychiatric genetics: the official publication of the International Society of Psychiatric Genetics. 2012;159b(1):30–7. pmid:22095641
          • 55. Veenstra-VanderWeele J, Xu T, Ruggiero AM, Anderson LR, Jones ST, Himle JA, et al. Functional studies and rare variant screening of SLC1A1/EAAC1 in males with obsessive-compulsive disorder. Psychiatric genetics. 2012;22(5):256–60. pmid:22617815
          • 56. Wendland JR, Moya PR, Timpano KR, Anavitarte AP, Kruse MR, Wheaton MG, et al. A haplotype containing quantitative trait loci for SLC1A1 gene expression and its association with obsessive-compulsive disorder. Archives of general psychiatry. 2009;66(4):408–16. pmid:19349310
          • 57. Gorinski N, Bijata M, Prasad S, Wirth A, Abdel Galil D, Zeug A, et al. Attenuated palmitoylation of serotonin receptor 5-HT1A affects receptor function and contributes to depression-like behaviors. Nat Commun. 2019;10(1):3924. pmid:31477731
          • 58. Clifton NE, Cameron D, Trent S, Sykes LH, Thomas KL, Hall J. Hippocampal Regulation of Postsynaptic Density Homer1 by Associative Learning. Neural plasticity. 2017;2017:5959182. pmid:29238619
          • 59. Urdánoz-Casado A, Sánchez-Ruiz de Gordoa J, Robles M, Acha B, Roldan M, Zelaya MV, et al. Gender-Dependent Deregulation of Linear and Circular RNA Variants of HOMER1 in the Entorhinal Cortex of Alzheimer’s Disease. Int J Mol Sci. 2021;22(17). pmid:34502114
          • 60. Bloom MA, Domm LV, Nalbandov AV, Bloom W. Medullary bone of laying chickens. American Journal of Anatomy. 1958;102(3):411–53. pmid:13617222
          • 61. Mueller WJ, Schraer R, Scharer H. Calcium metabolism and skeletal dynamics of laying pullets. The Journal of nutrition. 1964;84(1):20–6. pmid:14210016
          • 62. Johnsson M, Jonsson KB, Andersson L, Jensen P, Wright D. Genetic Regulation of Bone Metabolism in the Chicken: Similarities and Differences to Mammalian Systems. PLoS Genetics. 2015;11(5):e1005250. pmid:26023928
          • 63. Rybchyn MS, Brennan-Speranza TC, Mor D, Cheng Z, Chang W, Conigrave AD, et al. The mTORC2 Regulator Homer1 Modulates Protein Levels and Sub-Cellular Localization of the CaSR in Osteoblast-Lineage Cells. Int J Mol Sci. 2021;22(12). pmid:34204449
          • 64. Gajos-Michniewicz A, Piastowska AW, Russell JA, Ochedalski T. Follistatin as a potent regulator of bone metabolism. Biomarkers: biochemical indicators of exposure, response, and susceptibility to chemicals. 2010;15(7):563–74. pmid:20569048
          • 65. Koromila T, Dailiana Z, Samara S, Chassanidis C, Tzavara C, Patrinos GP, et al. Novel sequence variations in the CER1 gene are strongly associated with low bone mineral density and risk of osteoporotic fracture in postmenopausal women. Calcified tissue international. 2012;91(1):15–23. pmid:22543871
          • 66. Koromila T, Georgoulias P, Dailiana Z, Ntzani EE, Samara S, Chassanidis C, et al. CER1 gene variations associated with bone mineral density, bone markers, and early menopause in postmenopausal women. Human genomics. 2013;7(1):21. pmid:24138842
          • 67. Shin TH, Theodorou E, Holland C, Yamin R, Raggio CL, Giampietro PF, et al. TLE4 Is a Critical Mediator of Osteoblast and Runx2-Dependent Bone Development. Frontiers in cell and developmental biology. 2021;9:671029.
          • 68. van Karnebeek CD, Bonafé L, Wen XY, Tarailo-Graovac M, Balzano S, Royer-Bertrand B, et al. NANS-mediated synthesis of sialic acid is required for brain and skeletal development. Nat Genet. 2016;48(7):777–84. pmid:27213289
          • 69. Fullerton PT Jr., Monsivais D, Kommagani R, Matzuk MM. Follistatin is critical for mouse uterine receptivity and decidualization. Proc Natl Acad Sci U S A. 2017;114(24):E4772–e81.
          • 70. Liu Y, Fan J, Yan Y, Dang X, Zhao R, Xu Y, et al. JMY expression by Sertoli cells contributes to mediating spermatogenesis in mice. The FEBS journal. 2020;287(24):5478–97. pmid:32279424
          • 71. Sun SC, Sun QY, Kim NH. JMY is required for asymmetric division and cytokinesis in mouse oocytes. Molecular human reproduction. 2011;17(5):296–304. pmid:21266449
          • 72. Abramiuk M, Bębnowska D, Hrynkiewicz R, Niedźwiedzka-Rystwej P, Polak G, Kotarski J, et al. CLTA-4 expression is associated with the maintenance of chronic inflammation in endometriosis and infertility. Cells. 2021;10(3):487. pmid:33668701
          • 73. Iyer CC, Chugh D, Bobbili PJ, Iii AJB, Crum AE, Yi AF, et al. Follistatin-induced muscle hypertrophy in aged mice improves neuromuscular junction innervation and function. Neurobiology of aging. 2021;104:32–41. pmid:33964607
          • 74. Baker J, Liu J-P, Robertson EJ, Efstratiadis A. Role of insulin-like growth factors in embryonic and postnatal growth. Cell. 1993;75(1):73–82. pmid:8402902
          • 75. Mao J, Vanderlelie JJ, Perkins AV, Redman CW, Ahmadi KR, Rayman MP. Genetic polymorphisms that affect selenium status and response to selenium supplementation in United Kingdom pregnant women. The American journal of clinical nutrition. 2016;103(1):100–6. pmid:26675765
          • 76. Linke V, Overmyer KA, Miller IJ, Brademan DR, Hutchins PD, Trujillo EA, et al. A large-scale genome-lipid association map guides lipid identification. Nature metabolism. 2020;2(10):1149–62. pmid:32958938
          • 77. Catalán A, Merondun J, Knief U, Wolf JB. Chromatin accessibility, not 5mC methylation covaries with partial dosage compensation in crows. Plos Genetics. 2023;19(9):e1010901. pmid:37747941