Dosage Compensation throughout the Schistosoma mansoni Lifecycle: Specific Chromatin Landscape of the Z Chromosome

Abstract Differentiated sex chromosomes are accompanied by a difference in gene dose between X/Z-specific and autosomal genes. At the transcriptomic level, these sex-linked genes can lead to expression imbalance, or gene dosage can be compensated by epigenetic mechanisms and results into expression level equalization. Schistosoma mansoni has been previously described as a ZW species (i.e., female heterogamety, in opposition to XY male heterogametic species) with a partial dosage compensation, but underlying mechanisms are still unexplored. Here, we combine transcriptomic (RNA-Seq) and epigenetic data (ChIP-Seq against H3K4me3, H3K27me3, and H4K20me1 histone marks) in free larval cercariae and intravertebrate parasitic stages. For the first time, we describe differences in dosage compensation status in ZW females, depending on the parasitic status: free cercariae display global dosage compensation, whereas intravertebrate stages show a partial dosage compensation. We also highlight regional differences of gene expression along the Z chromosome in cercariae, but not in the intravertebrate stages. Finally, we feature a consistent permissive chromatin landscape of the Z chromosome in both sexes and stages. We argue that dosage compensation in schistosomes is characterized by chromatin remodeling mechanisms in the Z-specific region.


Introduction
Sex determination systems are very diverse and can involve genetic and/or epigenetic based mechanisms (Bachtrog et al. 2014). Genetic sex determination has been widely studied and often involves well-differentiated pairs of sex chromosomes: X and Y in male heterogametic systems, or Z and W in female-heterogametic systems. Morphological and gene content differences arise between the sex-specific Y/W and the shared X/Z chromosomes after their evolution from an ancestral pair of autosomes (Charlesworth 1991;Charlesworth et al. 2005). Successive events of recombination suppression (Rice 1987;Bergero and Charlesworth 2009) result in the accumulation of deleterious mutations and lead to the degeneration of the Y/W heterochromosome (Engelstadter 2008). Consequently, the sex carrying the degenerated Y/W harbors a certain number of genes in a single X/Z-linked copy. In the absence of mechanisms to buffer expression levels, such monosomy in a diploid genome is expected to induce detrimental effects on finely tuned gene networks (Veitia 2005). Mechanisms of gene expression regulation have evolved independently across eukaryotes to compensate for this gene dose imbalance and are grouped together under the name "dosage compensation." In the case of global dosage compensation, the overall expression level of monosomic genes on the heterochromosome (i.e., X/ Z-specific genes) is equal to the overall diploid expression level of autosomal genes (i.e., X/AA or Z/AA ratio of 1, "A" standing for the autosomal expression) (Vicoso and Bachtrog 2009). Conversely, the gene dose can be only partially compensated and the average expression of monosomic X/Zspecific genes will then be below the average expression of the diploid autosomes. This X/Z-specific gene expression value is typically not equal to half the dose of autosomal genes, as would be expected in the total absence of compensation, but instead reaches 60-80% of autosomal expression (i.e., X/AA or Z/AA ratio of 0.6-0.8). This is thought to reflect a combination of general buffering mechanisms (Stenberg et al. 2009) and/or the individual upregulation of dosage-sensitive sex-linked genes, that is, a "gene-by-gene" mechanism of compensation (Mank 2009). Global dosage compensation was initially thought to be not only 1) the rule but also 2) the preserve of XY systems. Both considerations are actually controversial, and many studies have challenged this canonical view. First, global compensation is not an absolute rule among XY species: initially, "chromosome-wide" compensation was described in both vertebrate and invertebrate XY model species, as well as a variety of non-model organisms  Darolti et al. 2019;and other vertebrates, reviewed in Graves 2016). Thereafter, the multiplicity of studies carried in a broad range of species actually highlighted a lack of global dosage compensation in vertebrates as well as in invertebrate (Xiong et al. 2010;Mank 2013;Hurst et al. 2015;Gu and Walters 2017 for review). Instead other strategies were observed, such as a geneby-gene mechanism with an increase of only some X-linked genes (Lin et al. 2012;Pessia et al. 2012;Albritton et al. 2014), the decreased expression of autosomal genes (Julien et al. 2012), or duplication-translocation of X-linked genes to autosomes (Hurst et al. 2015). Second, global compensation is not the preserve of XY systems: recent studies have shown that non-model ZW species also display total equalization of expression between Z-specific genes and autosomes (Lepidoptera: Walters and Hardcastle 2011;Kiuchi et al. 2014;Gu and Walters 2017;Huylmans et al. 2017; Artemia franciscana crustacean species: Huylmans et al. 2019). This is in contradiction to earlier studies in which only partial compensation was documented in many ZW female heterochromatic clades: birds (Itoh et al. 2007;Naurin et al. 2011;Wolf and Bryk 2011;Uebbing et al. 2015), arthropods (Harrison et al. 2012;Mahajan and Bachtrog 2015), snakes (Vicoso et al. 2013), flatfish Cynoglossus semilaevis (Chen et al. 2014), and metazoan parasites of the genus Schistosoma (Vicoso and Bachtrog 2011;Picard et al. 2018).
Schistosomes are blood flukes responsible for schistosomiasis, an infectious disease affecting more than 230 million people worldwide (Colley et al. 2014, for review). The model species Schistosoma mansoni has a complex life cycle, characterized by 1) clonal multiplication in a freshwater snail of the Biomphalaria genus, the intermediate host and 2) sexual reproduction in a definitive vertebrate host (i.e., a primate or rodent species). The parasite's eggs are released in freshwater via the feces. Free-living larvae (miracidia) hatch out and infect the mollusk intermediate host, transforming into sporocysts.
Sporocyst clonal multiplication inside the mollusk ultimately leads to thousands of infective cercariae, which are in turn released into fresh water. Definitive host penetration occurs through the epidermis and is followed by drastic morphological and physiological transformations: 1) within 2 hours, the free-living larvae become obligatory endoparasitic schistosomula; 2) after 2-5 weeks within the definitive host, these schistosomula develop from 150-lm juvenile sexually undifferentiated individuals into 1-cm adults. Schistosoma is the only genus displaying separate sexes among flatworms (Basch 1990;Combes 1991). Sex is genetically determined, with ZZ males and ZW females (Grossman et al. 1981), but no apparent phenotypic sexual dimorphism exists from the egg to the early stages of schistosomula. Phenotypic differences between males and females (i.e., sexual dimorphism) appear only in late intravertebrate stages (Loker and Brant 2006). In addition, male and female worms remain immature (i.e., sexually nonfunctional) until pairing. In particular, mated females undergo extensive morphological changes and develop their reproductive organs (Erasmus 1973;Neves et al. 2005).
Schistosomes, and particularly S. mansoni, have well-differentiated sex chromosomes: the W chromosome is mostly heterochromatic and carries many repetitive sequences (Portela et al. 2010;Lepesant et al. 2012;Protasio et al. 2012); 1,067 Z-specific genes have been identified thus far (Protasio et al. 2012;Picard et al. 2018). These Z-specific genes are consistently male-biased in expression, which was initially interpreted as absence of a global mechanism of dosage compensation (Vicoso and Bachtrog 2011). It was recently suggested that instead schistosomes have partially upregulated their Z chromosome in response to gene loss on the W (Picard et al. 2018), but that this upregulation occurred in both sexes (thereby partially balancing dosage in females, but maintaining the male-bias in expression). The conclusions of these studies were limited by several drawbacks. First, only older, vertebrate-infecting, stages were considered. Gene balance is thought to be most crucial in the earlier stages of development because of a higher number of functional interactions (Cutter and Ward 2005;Davis et al. 2005;Artieri et al. 2009;Mank and Ellegren 2009), and missing these early stages could lead to an underestimation of the extent to which compensation occurs. Similarly, male-biased expression can be due to genes with male-specific functions, which should largely be absent before the onset of sexual differentiation. Characterizing the sex-bias of Z-specific genes throughout development is therefore crucial to fully understand the dynamics of dosage compensation in this group. Second, fully distinguishing between local and global mechanisms of compensation requires an understanding of the molecular processes at work. In particular, well-studied chromosome-wide mechanisms of dosage compensation are characterized by global and extensive changes in the epigenetic landscape of the affected sex chromosome (Lucchesi et al. 2005;Lucchesi 2018). For instance, active chromatin marks are strongly enriched on the Drosophila male X, whereas the inactivated mammalian X chromosome of females is characterized by DNA methylation and widespread repressive chromatin marks (Lucchesi et al. 2005;Brockdorff and Turner 2015). On the chicken Z chromosome, a local and female-specific hyperacetylation of the fourth histone (H4K16Ac) has been described (Bisoni et al. 2005). Apart from this study, the detailed characterization of sex-linked histone chromatin marks is entirely lacking in ZW species, severely limiting our mechanistic understanding of dosage equalization in these groups.
Here, we systematically assess gene expression regulation along the Z chromosome of S. mansoni by combining transcriptomic (RNA-Seq) and epigenomic (ChIP-Seq) data. We focus on three different developmental stages: 1) cercariae (RNA-Seq and ChIP-Seq): these free larvae are the last fully sexually undifferentiated stage before host penetration; 2) schistosomula (RNA-Seq): the schistosomula were at an advanced stage of development but lacked phenotypic dimorphism (Picard et al. 2016); and 3) immature worms (RNA-Seq and ChIP-Seq): male and female parasite displaying sexual dimorphism, but sexually immature as they were not mated.

Publicly Available RNA-Seq and ChIP-Seq Reads
This study is mainly based on publicly available data that were previously generated in our laboratory. Transcriptomic data (i.e., total RNA isolation followed by single-read 50 nucleotides Illumina sequencing) for both sexes in cercariae, schistosomula (stage S#2), and immature worms were described in Picard et al. 2016 (raw reads are available under accession number SRP071285 on NCBI-SRA database). As for epigenetic data (i.e., chromatin immunoprecipitation, or ChIP, followed by single-read 50 nucleotides Illumina sequencing): H3K27me3 ChIP-Seq was analyzed in Picard et al. 2016 (male data available under accession number SRP071285); H3K27me3, H3K4me3, and H3K20me1 ChIP-Seq in females were reported in Roquis et al. 2016 andRoquis et al. 2018 (all data available on SRP035609). Accession numbers for each studied library are provided in table 1.

Newly Generated H3K4me3 ChIP-Seq Data in Males
Two biological replicates of male cercaria (2Â 10,000 individuals) and immature worms (2Â 20 individuals) were respectively obtained by monoclonal infection of Biomphalaria glabrata followed by unisex infection of Swiss OF1 mice (Picard et al. 2016 for details). Native chromatin immunoprecipitation assay was done according to Cosseau et al. (2009) using 4 ml of H3K4me3 antibody (Millipore, cat. number 04-745, lot number NG1680351). Further details are available at http://methdb.univ-perp.fr/epievo/; Last accessed June 27, 2019. ChIP library construction and sequencing were performed at the sequencing facilities of Montpellier GenomiX (MGX, France). Briefly, TruSeq ChIP sample preparation kit (Illumina Inc., USA) was used according to the manufacturer's recommendations on 30 ng of DNA per condition. DNAs were blunt ended and adenylated on 3 0 ends. Illumina's indexed adapters were ligated to both ends, and resulting ligated DNA were enriched by polymerase chain reaction (PCR). PCR products were separated by size using electrophoresis and 400 base pairs (bp) fragments were selected. The quantitative analysis of the DNA library was carried on Agilent High Sensitivity chip and qPCR (Applied Biosystems 7500, SYBR Green). Finally, the sequencing was performed on a HiSeq2500 in single-read 50-nt mode.

Gene Expression Analysis According to Genomic Location
Gene location was based on GTF annotations obtained at WormbaseParasite (schistosoma_mansoni.PRJEA36577.WB PS9.canonical_geneset.gtf, https://parasite.wormbase.org/index.html; Last accessed June 27, 2019). More precisely, genes located on the ZW sex chromosome linkage map were then attributed either to the Z-specific region or to pseudoautosomal regions (PARs) depending on their coordinates (Z-specific windows: Z1: 3,550,000-13,340,000; Z2: 13,860,000-19,650,000; Z3: 23,230,000-30,820,000, supplementary Data 3, Supplementary Material online) (Protasio et al. 2012). Normalized expression data combined to gene locations were then implemented into an R script for drawing gene expression patterns (supplementary Data 4, Supplementary Material online). This R script is available as supplementary Methods 1, Supplementary Material online.
Briefly, a Loess Normalization (R library Affy) was performed, taking into account the 12 libraries (cercariae, schistosomula, immature worms, for both sexes, and in duplicate). For each sex and stage, the averaged expression was then considered. Different thresholds for minimum of expression level were applied (RPKM > 0, RPKM > 1, and RPKM > 5), and strong sex-bias were considered or not. These strong sex-bias were defined by a fold change higher than 2, considering either the female-to-male ratio (female sex-bias), or the male-to-female ratio (male sex-bias). When comparing two conditions, the ratio between the medians of expression was calculated, and the level of significance was tested with Wilcoxon rank sum tests with continuity correction. In the figures, significance is showed by stars: *P value < 0.05, **P value < 0.001, and ***P value < 0.0001.

ChIP-Seq Raw Read Processing
ChIP-Seq data treatment was carried out under a local galaxy instance (Goecks et al. 2010). After quality check (Andrews 2010), neither quality filtering nor trimming was applied and all the reads were mapped to the S. mansoni reference genome (assembly version 5.2) (Protasio et al. 2012), using Bowtie2 (Langmead and Salzberg 2012). Mapping quality in Bowtie 2 is related to "uniqueness" of the read. SAM alignment files were converted into the bed format with pyicos (Althammer et al. 2011) and sorted with sortBed -i of the bedtools suite (Quinlan et al. 2011).

Comparative EpiChIP Analysis
Average histone modification profiles around transcriptional start site (TSS) of the genes were generated by doing a window analysis from À1,000 to þ5,000 base pairs relative to this TSS, using EpiChIP v0.9.7-e (Hebenstreit et al. 2011). As input, we used the 23 million, 6 million, and 19 million randomly sampled mapped reads that were generated after the alignment step for H3K27me3, H3K4me3, and H4K20me1, respectively. The average histone profiles were generated on the chromosome 1 and independently on the Z-specific region and the PAR of the ZW sex chromosomes. For this purpose, we used the GTF annotation file generated previously (Picard et al. 2016) (available under id "Schistosoma mansoni sex-specific transcriptome" at http://ihpe.univ-perp.fr/accesaux-donnees/; Last accessed June 27, 2019) and selected 6,225 transcripts on chromosome 1, 2,421 transcripts (¼ 2421) on the Z-specific region, and 3,376 on the PAR. The average H3K4me3, H3K27me3, and input profiles (i.e., control without antibody) were generated for the two male biological replicates and the three female biological replicates (Supplementary Data 6 and 7, Supplementary Material online). The average H4K20me1, and input profiles were generated for the three female biological replicates (Supplementary Data 8, Supplementary Material online). Each average profile was normalized with its respective input average profile in order to be able to compare the different regions, whatever the number of genes in the considered -RNA-Seq and ChIP-Seq data were previously generated in our laboratory and were published in two different studies: SRP071285 accession number (SRA-NCBI database) for RNA-Seq and male H3K27me3 ChIP-Seq duplicates (Picard et al. 2016); and SRP035609 accession number (SRA-NCBI database) for female ChIP-Seq triplicates (Roquis et al. 2016). H3K4me3 ChIP-Seq analysis on males was never presented before and was deposited under SRP071285 accession number (SRA-NCBI database).
region. The distribution of histone enrichment around the transcription start site was compared according to the stage and the genomic location using Kolmogorov-Smirnov two sample tests.
Code Availability R script for expression analysis is provided in supplementary Methods 1, Supplementary Material online.

Different Levels of Gene Expression Equalization throughout Development
Gene expression was analyzed by using published RNA-Seq data from male and female S. mansoni of three developmental stages, characterized by different parasitic status and levels of sexual differentiation (Picard et al. 2016): cercariae, schistosomula, and immature worms (Fig. 1A ). Other stages such as miracidia or eggs were excluded because it is not possible to distinguish males and females at this phase. To assess dosage compensation, we compared both the Z-to-autosome ratio within each sex (Z:AA in females or ZZ:AA in males), and the female-to-male ratio between sexes (Z[F:M]/A[F:M]) (table 2). Such comparisons support global dosage compensation if the female Z:AA ratio and the Z(F:M)/A(F:M) ratio are equal to one; lower ratios support a lack of global dosage compensation. Different minimum expression thresholds were applied, and strongly sex-biased genes (>2-fold difference between the sexes) were excluded, as the presence of genes with sex-specific functions can also lead to male-biased expression of sex chromosomes, even in the presence of global dosage compensation (Huylmans et al. 2017). Gene expression ratios and their level of significance are detailed in table 2 for the three stages and six different methodological conditions. In order to minimize noise but keep a reasonable sample size, we focused on results obtained with a minimum expression threshold of RPKM > 1 and excluding genes with a sex-bias greater than 2-fold (Sample sizes are reported in cercariae are free larvae without any sexual dimorphism (1); schistosomula represent an intravertebrate stage, with ongoing sexual differentiation but no phenotypic sexual dimorphism (2); and immature worms are sexually differentiated but sexually non-functional as they did not mate (3). Expression of autosomal genes ("A," dark shade), pseudoautosomal genes ("PAR," dark shade), and Z-specific genes ("Z," light shade) is represented considering femaleto-male ratio ("F:M") (B), or independently within females (C) and males (D). Only genes with expression RPKM > 1, and a female-to-male fold change lower than 2 (sex-bias filtering) are taken into account (n ¼ 3,741 in cercariae; n ¼ 5,636 in schistosomula; n ¼ 5,657 in immature worms). The Z-to-autosome expression ratio for each sex (Z:AA for female, and ZZ:AA for male), and the corresponding female-to-male ratio (Z[F:M]/A[F:M]) are detailed in table 2 ("RPKM > 1 -sex-bias filtered"). Asterisks show the level of significance for each of these comparisons (Wilcoxon test): *P value < 0.05, **P value < 0.001, and ***P value < 0.0001.
Supplementary Table 1, Supplementary Material online), but qualitative patterns generally hold for other filtering procedures.
As observed before (Vicoso and Bachtrog 2011;Picard et al. 2018), the female-to-male ratio of expression was always significantly lower on the Z than on the autosomes, and this seemed to be due to both reduced Z:AA expression in females (Z:AA from 0.75 to 1.05) and an excess of ZZ:AA expression in males (ZZ:AA from 1.07 to 1.75) (Picard et al. 2018). However, there were important differences between the developmental stages. First, the Z(F:M)/A(F:M) ratio ranged from 0.60 to 0.90 (table 2C). Although part of this variation was driven by the different filtering procedures, in every case the value was closer to 1 for cercariae, consistent with more extensive dosage compensation in early development. For instance, for Z-specific genes with RPKM > 1 and no strong sex-bias, Z(F:M)/A(F:M) was equal to 0.81, 0.70, and 0.64 in cercariae, schistosomula, and immature worms, respectively, with significant differences between the free larval stage cercariae and the intravertebrate stages (P value < 0.0001, fig. 1B).
This difference in the extent of dosage equalization was also apparent within one sex. In the females, discrepancies were observed depending on the parasitic status and the filtering method: for RPKM > 0, after sex-bias filtering, the Z-toautosome ratio (Z:AA) ranged from 0.75 in female schistosomula (lower than 1, consistent with partial dosage compensation) to 1.05 in female cercariae (equalized ratio, typical of dosage global compensation) and was not always significant from 1 (table 2A). This seemed to be largely driven by the dichotomy between cercariae and intravertebrate stages ( fig. 1C): when considering genes with RPKM values above 1, and after filtering for sex-biased genes, Z-specific genes and autosomes displayed an equalized expression in female cercariae (Z:AA ¼ 1.00 for RPKM > 1, Z:AA ¼ 1.03 for RPKM > 5), but not in female schistosomula or immature worms (Z:AA ¼ 0.76 and 0.87 for RPKM > 1, Z:AA ¼ 0.90 and 0.80 for RPKM > 5). Although the specific numbers varied, cercariae had the highest female Z:AA expression ratio for all but one filtering procedures (RPKM > 0 and no removal of strongly sex-biased genes), generally supporting the more extensive upregulation of expression of Z-linked genes in females at this stage. Within males, the Z-to-autosome ratio (ZZ:AA) was higher in all studied conditions (ranging from 1.07 to 1.75) although the level of significance varied depending on the filtering procedure (table 2B). For instance, considering genes with RPKM > 1 and no strong sex-bias, the ZZ:AA ratios are 1.07, 1.12, and 1.27 for cercariae, schistosomula, and immature worms, respectively ( fig. 1D).
In summary, this first transcriptomic analysis of Z-specific and autosomal genes throughout schistosome development showed that 1) the overexpression of the male Z-specific genes was found consistently in the three stages; 2) the female Z-to-autosome ratio was around 1 in cercariae, as expected for global dosage compensation; and 3) the female Z-to-autosome ratio was around 0.8 in the intravertebrate stages, as expected when partial dosage compensation occurs.

Regional Variation of Dosage Compensation along the Z-Specific Region
In some species which display partial dosage compensation, the female-to-male ratio has been shown to vary along the Z chromosome, with some regions fully compensated, whereas others are not compensated at all (e.g., Gallus gallus  ). To investigate potential regional variation in S. mansoni, we investigated local patterns of gene expression by using sliding windows of 50 genes along the Z chromosome ( fig. 2). The Z-specific part, which was previously described as discontinuous in the version 5.2 of the genome (Protasio et al. 2012), was considered here as three Z-specific regions named Z1, Z2, and Z3. They were respectively defined by their coordinates, as follow: Z1 from 3,550 to 13,340 kb; Z2 from 13,860 to 19,650 kb; and Z3 from 23,230 to 30,820 kb (Protasio et al. 2012). As expected, the sliding window analysis revealed a female-to-male ratio oscillating around 1 in the PAR for all stages (log 2[F:M] close to 0; fig. 2, right panel). In the three Zspecific regions, the female-to-male expression ratio oscillated around 0.7 in schistosomula and 0.6 immature worms, consistent with homogeneous but partial compensation ( fig. 2D and F; supplementary figs. 1 and 2, Supplementary Material online). However, in cercariae, the female-to-male expression ratio was closer to 1. When looking at each Z-specific region individually, the female-to-male expression ratio was significantly lower in the Z1 and Z3 regions compared with the PAR in all three developmental stages ( fig. 2 fig. 2A). In schistosomula and immature worms, the same Z2 was significantly more male biased in expression than the PAR (Z2[F:M]/PAR[F:M] ¼ 0.71 and 0.62, respectively; fig. 2C and E).
Consistent with this, male and female expression levels overlapped in the Z2 region of cercariae; whereas in the Z1, female expression was lower than male expression, consistent with less complete compensation ( fig. 2B). In the Z3 region, and still in cercariae, we found a narrow window of overlapping male and female expression at the beginning of the region, whereas the remaining part displayed lower S. for Z1 < Z3 < Z2), which could be driven by these smaller-scale regional differences ( fig. 2B). In the intravertebrate stages, the three Z-specific regions differed clearly from the pseudoautosomal ones by presenting consistently reduced expression in females ( fig. 3D and F). The female and male gene expression ratios for each Z-specific region are shown in Supplementary Table 2 In summary, this transcriptomic analysis along the Z chromosome highlighted 1) local variation of the dosage compensation across Z-specific regions in cercariae and 2) partial and consistent dosage compensation all along the three Z-specific regions in schistosomula and immature worms.

Z-Specific Chromatin Landscape in Cercariae and Immature Worms to Elucidate Shared or Sex-Specific Mechanisms
Previous studies on dosage compensation in model organisms indicate that chromatin structure plays a key role in the regulation of this evolutionary mechanism (Lucchesi et al. 2005;Lucchesi 2018). As the main transcriptomic differences were observed between the free stage cercariae and the intravertebrate stages, we focused our epigenetic study on cercariae and immature worms. We analyzed immunoprecipitation assays with antibodies targeting three histone marks: 1) H3K4me3, associated with active promoter and transcription start site, and depleted in the heterochromatic inactivated X chromosome of mammals (O'Neill et al. 2008;Marks et al. 2009); 2) H3K27me3, a repressive mark associated with polycomb, poised transcription, and enriched in the inactivated X chromosome of mammals (Lucchesi et al. 2005); and 3) H4K20me1, associated with nonpermissive chromatin in Caenorhabditis dosage compensation (Vielle et al. 2012;Kramer 2015Kramer , 2016. We compared the enrichment of these histone marks in a sex-specific manner, on Z-specific genes, pseudoautosomes, and autosomes (using chromosome 1 as a proxy for the autosomal chromatin landscape, as its transcriptomic pattern appeared to be representative of autosomes, see Supplementary Figure 5, Supplementary Material online). For each of these genomic locations, the average enrichment profile of all annotated genes was performed around the TSS: from 1,000 base pairs (bp) upstream to 5,000-bp downstream.
We observed that the average profile of the three studied marks was specific for each genomic location in both cercariae and immature worms (Figure 3 and supplementary figs. 6-8, Supplementary Material online). Z-specific genes were enriched for H3K4me3 in both stages and sexes, especially between the TSS and position þ2,000 bp ( fig. 3, upper panel). H3K27me3 was depleted from Z-specific genes upstream of FIG. 3.-Average H3K4me3 and H3K27me3 enrichment profile according to sex, developmental stage, and genomic location. x axis represents the position in base pairs (bp) relative to the TSS of the genes (position 0). y axis represents the normalized average enrichment of reads obtained after a Chromatin Immunoprecipitation targeting the "permissive" mark H3K4me3 (in green) and the "nonpermissive" H3K27me3 (in red), in female cercariae (A) and immature worms (B), or male cercariae (C) and immature worms (D). The EpiChIP enrichment has been calculated around the TSS for chromosome 1 as proxy for autosomes (Chr1, dark shade), for the PAR of sex chromosomes (PAR, medium shade), and for chromosome Z-specific region of sex chromosomes (Z, light shade). For each of these genomic locations, we show the average result of the profiles obtained for each coding sequence. Each profile has been normalized with the same average enrichment of reads obtained after a Chromatin Immunoprecipitation without antibody. The experiment was performed in duplicates in males and triplicate in females. EpiChIP profiles showing standard error at each position are shown in Supplementary Figures 6 and 7, Supplementary Material online. The percentage of maximum difference between genomic regions is shown in Supplementary Tables 3 and 4, Supplementary Material online: all differences are statistically significant (P value < 0.001, Kolmogorov-Smirnov two sample tests).
the TSS and along the transcription unit in both stages and sexes ( fig. 3, lower panel). Finally, H4K20me1, for which only female data was available, displayed a stronger enrichment on chromosome 1 and on the PAR than on Z-specific genes in female cercariae (percentage of maximum difference: 63.3%, supplementary fig. 8 and supplementary table 5, Supplementary Material online). However, this difference diminished between stages, with an increased enrichment on the Z-specific region in immature females bringing it closer to the autosomal level (percentage of maximum difference: 30%, supplementary fig. 8 and supplementary table 5, Supplementary Material online). Differences between chromosomal regions within a developmental stage and sex were tested for the three analyzed histone marks: the most significant difference was consistently observed either between Z-specific genes and chromosome 1 or between Zspecific and pseudoautosomal genes (Kolmogorov-Smirnov two sample tests, supplementary tables 3-5, Supplementary Material online). This supports the idea that the chromatin landscape is more similar between the autosomes and the PAR of S. mansoni, whereas Z-specific genes display a singular epigenetic landscape. The depletion of the two "repressive" marks (H3K27me3 and H4K20me1) and the enrichment of the "permissive" H3K4me3 mark suggest a chromatin structure prone to enhanced gene expression specifically in this Zspecific region.
In summary, our epigenetic analysis highlighted the specialized chromatin landscape of the Z-specific region, which appeared to be more permissive than that of the autosomal regions in both sexes and stages. This modified chromatin state of the Z-specific region likely promotes enhanced expression compared with the autosomes, independent of the stage and sex, consistent with the first part of our study. The analysis of H4K20me1 in females further showed a differential enrichment of this histone mark on Z-specific genes between cercariae and adults, which may be related to the differences that we observe in the extent of female dosage compensation between the two stages.

Discussion
We provide here the first combined transcriptomic and epigenetic analysis of Z chromosome dosage compensation through three developmental stages of S. mansoni parasite. This work highlights three important aspects of gene expression regulation on sex chromosomes: 1) A strong upregulation of female expression of Z-linked genes in the free larval stage cercariae, consistent with complete dosage compensation, but not in intravertebrate stages schistosomula and immature females that display partial dosage compensation only. 2) Local variations of the female-to-male expression ratio along the Z chromosome. 3) Differences in chromatin structures between Z-specific regions and autosomes which may support the enhanced expression of both male and female Z chromosomes.

Change in Dosage Compensation Status following Host Penetration
Our RNA-Seq analysis detects a global hypertranscription of Zspecific genes in both sexes of S. mansoni, consistent with the intermediate stage of compensation previously reported in schistosomula and adult worms (Vicoso and Bachtrog 2011;Picard et al. 2018). In both studies, the female upregulation of expression only partially resolved the imbalance between the Z chromosome and the autosomes. However, we show here that this partial compensation does not apply to earlier developmental stages of S. mansoni: in females, enhanced transcriptional levels of the Z-specific genes results in a Z-to-A ratio of around 1 in cercariae (compared with $0.8 in the intravertebrate stages).
The dosage compensation status has been shown to be tissue-and development-specific within a same species (Mank and Ellegren 2009;Lott et al. 2011;Deng et al. 2014;Huylmans et al. 2017). In particular, germ/stem cells might need to overcome a special challenge regarding dosage compensation, as 1) genome wide reprograming is expected to occur and erase epigenetic marks responsible for maintaining the chromatin state necessary for compensation (Sangrithi and Turner 2018) and/or 2) dosage compensation mechanisms may act as developmental regulators by targeting autosomal genes involved in patterning and morphogenesis (Valsecchi et al. 2018). In schistosomes, larval stages display particular features regarding their content in such cells: in the intermediate host, sporocysts consist of totipotent stem cells undergoing intense clonal multiplication (Cort et al. 1954); in cercariae, embryonic stem cell specific combination of histone marks have been described (Wang et al. 2013;Roquis et al. 2015). This chromatin signature disappears after host penetration, signifying important changes in gene expression regulation (Roquis et al. 2015). Therefore, a specific cell type content could explain the shift in the dosage compensation status observed between the larval stage and the intravertebrate phase. Further studies are needed to disentangle the role of those cells in the observed pattern, and the examination of dosage compensation status during earlier larval stages certainly deserves further attention. Technological advances such a single cell/single individual RNA sequencing should in the future allow for such analyses.
Local Variation of Male-Bias in Cercariae: Different Evolutionary Stages of the Dosage Compensation Mechanisms?
Overexpression of Z-specific genes compared with the level of autosomal gene expression has been recently described in S. mansoni intravertebrate stages (Picard et al. 2018). Here, we further show that this overexpression also occurs in the undifferentiated and free-living stage cercariae and is therefore a consistent feature maintained throughout the parasite life cycle. This uncommon feature has been described in adult floor beetles (whole body sample) which is a XY male heterogametic system (Prince et al. 2010). In this species, the male hemizygous X chromosome is fully compensated, and females display XX/AA ratio superior to 1. It has been also recently reported in the young XY plant system Silene latifolia (Muyle et al. 2018). Such pattern was theorized by Gu and Walters as the "type IV" of sex chromosome dosage compensation, with complete equalization of the X/Z-toautosome ratio in the heterogametic sex, but no dosage balance between sexes as the homogametic sex exhibit Z/X hypertranscription relative to autosomes (Gu and Walters 2017). These observations can be interpreted as an example of the first step of the Ohno's hypothesis, suggesting that the establishment of dosage compensation process evolved in a two-step process (Ohno 1967;Charlesworth 1996;Vicoso and Bachtrog 2009;Bachtrog et al. 2011;Mank 2013): 1) First, Z-or X-linked genes became upregulated in both sexes. This would restore expression of the single X or Z chromosome in the heterogametic sex to the diploid level that existed before degradation of the Y or W, but would result in too much gene expression in the homogametic sex. 2) Therefore, secondarily, a downregulation of Z or X-linked gene expression would evolve in the homogametic sex to restore the correct Z-or X-to-autosome ratio. But many studies recently challenge this theory. For instance, there is no convincing evidence that the X chromosome is globally upregulated in both mammals (Julien et al. 2012), and C. elegans (Albritton et al. 2014). In C. elegans, the X chromosomes in the XX homogametic individuals are downregulated, leading to dosage balance between the sexes but not dosage compensation in X0 males. In placental mammals, no global X upregulation was observed, unlike in marsupials (Julien et al. 2012;Whitworth and Pask 2016). Thus, even though complete X chromosome inactivation has traditionally been interpreted as the second step of the process (Brockdorff and Turner 2015), the relevance of the Ohno's hypothesis to the evolution of placental mammal dosage compensation is in fact controversial (e.g., Nguyen and Disteche 2006;Vicoso and Bachtrog 2009;Lin et al. 2012;Mank 2013;Pessia et al. 2014;Gu and Walters 2017).
Our results suggest that S. mansoni may have followed an "Ohno-like" evolutionary trajectory for dosage compensation, which allowed the Z-specific region to enhance its expression to reach that of autosomal genes in females, but under which males somehow avoid the necessity of complete countercompensation of this Z chromosome hypertranscription. This raises the question of what the current status for the evolution of dosage compensation is in S. mansoni: an intermediate or a stable state?
Regional variation of the Z dosage compensation has been described in birds, where some dosage-compensated genes are concentrated in a region of the short arm of the Z chromosome, near the male hypermethylated (MHM) locus (Melamed and Arnold 2007;Wright et al. 2015). Full dosage compensation has also been described in a restricted chromosomal region in pseudomale testes of Cy. semilaevis, whereas the rest of the Z chromosome is partially compensated ). It has been proposed that such local variations in gene regulation along the sex chromosome could be based on regional age following sex-linkage, including in XY species such as stickleback (Schultheiß et al. 2015). In S. mansoni, two evolutionary strata of different ages were recently described (supplementary fig. 9, Supplementary Material online) (Picard et al. 2018). If partial upregulation of Z-specific genes in both sexes represents an intermediate step in the evolution of dosage compensation, we may expect that the older stratum will be closer to balanced dosage between males and females than the younger one. In the previous study and in ours, no difference in gene expression could be detected between these strata in intravertebrate stages. In cercariae, however, the section of the Z that is closest to complete compensation (region Z2) is indeed part of the older Z-specific evolutionary stratum that is shared between African and Asian schistosomes (supplementary fig. 9, Supplementary Material online, Picard et al. 2018). It is therefore likely that dosage compensation is still evolving in S. mansoni, and that the Z-specific region which displays equalized gene expression between the sexes actually represents the final stage of dosage compensation evolution in this species.

The Specific Chromatin Features of the Z-Specific Regions May Account for the Hypertranscription of Z-Specific Genes in Both Sexes
Various mechanisms have evolved to regulate the gene dosage at the functional level. In fully compensated organisms, these mechanisms are all based on the modulation of chromatin accessibility of the sex-specific regions (Lucchesi et al. 2005;Ercan 2015;Lucchesi 2018). In mammals, transcriptional regulation is achieved by the progressive depletion of histone active marks concomitant to the enrichment of histone repressive mark H3K27me3 through the action of the Xinactive specific transcript (Xist) (Brockdorff and Turner 2015). In C. elegans (XX/XO system), the hermaphrodite two X chromosomes display halved transcription level to match X expression to that of males (XO). Reduced transcription is allowed by a complex of proteins called the dosage compensation complex resulting in a depletion of histone active marks and enrichment of histone repressive marks (Lau and Csankovszki 2015). In Drosophila melanogaster (XX/XY system), the transcription of the single male X chromosome is doubled by an overall increase of the chromatin accessibility by the Male-Specific Lethal complex (Lucchesi et al. 2005). More recent work performed on the pea aphid has also shown an enhanced chromatin accessibility of the X chromosome of males which may account for dose correction of those X-linked genes in this species (Richard et al. 2017).
We present here an overview of the chromatin structure by ChIP-Seq analysis targeting three modified histones in different developmental stages of S. mansoni and highlight different chromatin patterns between autosomal and Z-linked genes. We found that H3K27me3 and H4K20me1 were depleted in the Z-specific region relative to the PAR of Z chromosome and to autosomal regions. Reversely, we found that H3K4me3 was enriched in the Z-specific region of the Z chromosome relative to the PAR and to autosomal regions. The role of H4K20me1 in gene regulation has remained a mystery because of its contribution to both gene activation and gene repression in different contexts (Beck et al. 2012). However, in the context of chromosome-wide regulatory mechanism for dosage compensation, this mark has been shown to be enriched in the inactive chromosome of female mammals (Kelsey et al. 2015) and enriched in both hermaphrodite X chromosomes in Caenorhabditis, resulting in reduction of transcription level (Bian et al. 2017;Kramer 2015Kramer , 2016. Here, changes in H4K20me1 between female cercariae and adults are specifically observed for the Z-specific region, arguing in favor of a role of this mark in dosage compensation status switch between the stages. H3K27me3 is a clear repressive mark enriched in heterochromatic X-inactivated chromosome of mammals (Kelsey et al. 2015) and H3K4me3 was reported to be depleted during X-inactivation in female embryonic stem cells (O'Neill et al. 2008). Given the depletion for both repressive marks concomitant to enrichment of the active H3K4me3 mark in the Z-specific region, we suggest that a global regulation of the chromatin accessibility of the Z-linked regions occurs in order to overexpress the Z-linked genes and compensate for gene dose defect in the single Z-specific region of females. This could be further addressed using chromatin accessibility assay which could evidence an open chromatin state such as those performed in the pea aphids which support their enhanced X overexpression (Richard et al. 2017).
Local variation in dosage compensation is also supported by chromatin based event such as those described in C. semilaevis where an increase in cytosine methylation density occurs in the compensated region . In Gallus gallus, the implication of the MHM noncoding RNA and subsequent enrichment of H4K16ac around the MHM locus in females allow a full compensation of this region (Melamed and Arnold 2007). Female-specific non-coding RNAs have been shown to be expressed in schistosome larval stages (Lepesant et al. 2012) and their implication in changes in chromatin states, regulation of gene expression and, more specifically, in dosage compensation mechanisms, certainly deserve further attention.

Conclusion
Until recently, dosage compensation in ZW femaleheterogametic species was thought to be partial and to occur at a gene-specific level. Recent studies in non-model organisms have challenged this canonical view. In line with them, our study brings a new insight by showing developmental changes in dosage compensation status in female schistosomes. From global compensation in undifferentiated free larvae, to partial compensation after host penetration and the onset of sexual differentiation. Despite this developmental variation, our study highlights a robust overexpression of the Z chromosome throughout S. mansoni life cycle, independently of the sex. We show how this might be mediated by an enhanced chromatin accessibility of the Z-specific regions, giving a first insight into S. mansoni chromatin pattern in relation to dosage compensation. Our epigenetic study paves the way toward the construction of an evolutionary chromatin landscape of the parasite's dosage compensation. Investigating more combined histone modifications and non-coding RNAs appears to be the next step to understand both developmental changes and finer variations in gene expression along the Z chromosome.

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.

Ethic Statement
All experiments were carried out following the national ethical standards established in the writ of February 1, 2013 (NOR: AGRG1238753A), which set the conditions for approval, planning, and operation of establishments, breeders, and suppliers of animals used for scientific purposes and controls. The French Ministè re de l'Agriculture et de la Pêche and the French Ministè re de l'Education Nationale de la Recherche et de la Technologie provided permit A66040 to the laboratory for animal experiments and certificate to the experimenters (authorization 007083, decree 87-848).