Sex- specific expression and DNA methylation in a species with extreme sexual dimorphism and paternal genome elimination

Phenotypic differences between sexes are often mediated by differential expression and alternative splicing of genes. However, the mechanisms that regulate these expression and splicing patterns remain poorly understood. The mealybug, Planococcus

stages (Grath & Parsch, 2016). However, the mechanisms that regulate these sex-specific expression patterns are often poorly understood.
DNA methylation is a well-characterized epigenetic modification that could facilitate such variation in expression (Grath & Parsch, 2016). DNA methylation is found throughout the genome of many organisms (Suzuki & Bird, 2008) and occurs most frequently at 5′-CG-3′ dinucleotides, known as CpG dinucleotides (Bird, 1986). In mammalian somatic tissue, 70%-80% of all CpG sites are methylated (Feng et al., 2010) and methylation at promoter regions can suppress gene transcription, leading to stable gene silencing (Bird, 2002). This is implicated in the regulation of sex-specific and sex-biased gene expression (examples include Hall et al., 2014;Maschietto et al., 2017). In contrast, DNA methylation levels in arthropods are generally much sparser and vary across taxa (Thomas et al., 2020). In most holometabolous insects, DNA methylation is almost exclusively restricted to exons in a small subset of transcribed genes (Zemach et al., 2010). Overall the highest levels of global DNA methylation are found in hemimetabolous insects (e.g. 14% in Blattodea, Bewick et al., 2016), while methylation is much more limited in holometabolous species (Lewis et al., 2020;Provataris et al., 2018). In insects, the role of DNA methylation in the regulation of gene expression remains inconclusive. However, studies show that DNA methylation is generally associated with elevated, stable gene expression (Bonasio et al., 2012;Foret et al., 2009;Glastad et al., 2016;Wang et al., 2013).
It is also worth noting there is mounting evidence for the role of DNA methylation in alternate genome functions within arthropods. For example, Lewis et al. (2020) have recently shown that exon methylation across a number of arthropod species appears to be associated with nucleosome positioning. Additionally, recent evidence in honeybees suggests a possible role for DNA methylation in the maintenance of gene duplications (Dyson & Goodisman, 2020;Yagound et al., 2020);;; have found DNA methylation profiles tends to be associated with the underlying genotype of an individual. These studies highlight the possibility of multiple roles of DNA methylation within arthropod species and may explain why there is currently no conclusive evidence for whether DNA methylation can directly regulate gene expression in insects.
Despite some evidence suggesting a relationship between DNA methylation and gene expression, few insect studies have directly explored sex-specific DNA methylation patterns and their association with sex-specific gene expression. In the jewel wasp, Nasonia vitripennis, 75% of expressed genes show sex-biased expression; however, DNA methylation patterns between the sexes are similar and do not explain gene expression patterns (Wang et al., 2015). In contrast, a study in the peach aphid, Myzus persicae, in which 19% of genes exhibit sex-specific expression biases, reveals a positive correlation between changes in gene expression and DNA methylation, where higher sex-specific gene expression is associated with higher sex-specific methylation, particularly for genes located on the sex chromosomes (Mathers et al., 2019). However Glastad et al. (2016) find few DNA methylation differences between sexes in a termite species (Zootermopsis nevadensis) compared to the number of differences between castes. Thus, the role of sex-specific patterns of methylation in regulating sex-biased gene expression in insects remains unclear.
The citrus mealybug, Planococcus citri (Hemiptera: Pseudococcidae), is uniquely suited for studying the functional role of DNA methylation in sex-specific gene expression. In addition to extreme sexual dimorphism ( Figure 1), P. citri also has an unusual reproductive strategy, known as paternal genome elimination (PGE).
PGE is a genomic imprinting phenomenon found in thousands of insect species (Burt & Trivers, 2006) that involves the silencing and elimination of an entire haploid genome in a parent-of-origin specific manner. Under PGE, both sexes develop from fertilized eggs and initially possess a diploid euchromatic chromosome complement. However, males subsequently eliminate paternally inherited chromosomes during spermatogenesis and only transmit maternally inherited chromosomes to their offspring (Brown & Nelson-Rees, 1961). Furthermore, in P. citri males, paternally inherited chromosomes are heterochromatinized in early development (Bongiorni et al., 2001;Brown & Nur, 1964), and thus, gene expression shows a maternal bias (de la Filia et al., 2020). Females, on the other hand, do not undergo the process of PGE and both maternally and paternally derived chromosomes remain euchromatic throughout development (Brown & Nur, 1964). Due to the haploidization of males, PGE is often referred to as a 'pseudohaplodiploid' system. Furthermore, we have previously shown P. citri females have a unique pattern of whole-genome DNA methylation that differs from that found in other arthropods (Lewis et al., 2020). While most arthropods have depleted levels of transposable element and promoter methylation, P. citri has independently evolved both (Lewis et al., 2020). Interestingly, and similar to patterns shown in mammals, genes with low expression in P. citri have significantly higher promoter methylation than highly expressed genes (Lewis et al., 2020).
It is also suggested that DNA methylation may have a role in the recognition and silencing of paternally derived chromosomes in males in the process of PGE (Bongiorni et al., 1999;Buglia et al., 1999).
Supporting the idea that DNA methylation may be involved in sexual dimorphism and PGE in mealybugs and other scale insects, two recent studies have identified sex-biased expression of the DNA methyltransferase DNMT1 in adult Phenacoccus solenopsis (Omar et al., 2020) and Ericerus pela (Yang et al., 2015), with females showing considerably higher expression compared to males in both species.
In order to identify sex-specific patterns of gene expression and clarify the role of DNA methylation in this process, we analyse both male and female P. citri methylomes and transcriptomes. This is the first genome-wide analysis of sex-specific gene expression and DNA methylation in scale insects. Using RNA-seq and whole-genome bisulphite sequencing (WGBS), we find clear differences in gene expression and methylation profiles between the sexes. However, we find no relationship between differentially expressed genes and differentially methylated genes, indicating that cis-acting DNA methylation does not drive sex-specific gene expression in adult P. citri.

| Insect husbandry
Mealybug cultures used for this study were kept on sprouting potatoes in sealed plastic bottles at 25°C and 70% relative humidity. Under these conditions, P. citri has a generation time (time from oviposition until sexual maturity) of approximately 30 days.
Experimental isofemale lines were reared in the laboratory under a sib-mating regime: in each generation, one mated female is taken per culture and transferred to a new container to give rise to the next generation. The P. citri line used (WYE 3-2) was obtained from the pest control company, WyeBugs in 2011, and had undergone 32 generations of sib-mating prior to this experiment. This high degree of inbreeding allows for precise mapping of whole-genome Bisulfiteseq (WGBS) reads reducing mis-mapping caused by SNP variation.
It also means we avoid contrasting methylation profiles caused by differences in the underlying genotype of individuals (epialleles).
We isolated virgin females after they became distinguishable from males (3rd-4th instar) and kept them in separate containers until sexual maturity (>35-days old). Males were isolated at the pupal stage and kept in separate containers until eclosion (~27 days).
Insects were stored at −80°C until DNA and RNA extraction.

| Differential expression and alternative splicing
Raw RNA-seq reads for each sample were trimmed for low-quality bases and adapters using Fastp for paired-end reads (Chen et al., 2018).
Fastp was used as it allows removal of poly-G tails from NovaSeq reads. We quantified gene-level expression for each sample using rsem v1.2.31 (Li & Dewey, 2011) with star v2.5.2a (Dobin et al., 2016) based on the P. citri reference genome and annotation (mealybug.org, version v0). Additional information on the reference genome can be found in ing for multiple testing using the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995).
Illumina reads were first aligned to the converted unmethylated lambda E. coli control DNA sequence using bismark v0.19.0 (Krueger & Andrews, 2011) to estimate the error rate of the C-T conversion. bismark v0.19.0 and bowtie2 were then used to align reads to the reference genome using standard parameters. The weighted methylation level of each genomic feature (P. citri v0 annotation, mealybug.org) was calculated as in Schultz et al. (2012). Briefly, this method accounts for the CpG density of a region by calculating the sum of all cytosine calls for every CpG position in a region (promoter/exon/gene etc.) divided by the total cytosine and thymine calls in the same region.
For differential methylation analysis between sexes, coverage outliers (above the 99.9% percentile) and bases covered by <10 reads were removed. Each CpG per sample was subjected to a binomial test to determine the methylation state, where the lambda conversion rate was used as the probability of success. Only CpGs which were determined as methylated in at least one sample were then tested via a logistic regression model, implemented using methylkit v1.10.0 (Akalin et al., 2012), for differential methylation between the sexes. p-values were corrected for multiple testing using the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995).
CpGs were considered differentially methylated if they had a q-value <0.01 and a minimum methylation difference of 15%.
Promoter and exon regions were classed as differentially methylated if they contained at least three significant differentially methylated CpG sites and had a weighted methylation difference >15% across the entire region. Significant overlap of genes with promoter and exon differential methylation was determined using the hypergeometric test and visualized using the upsetr package v1.4.0 (Lex et al., 2016).

| Relationship of gene expression and DNA methylation
The relationship of promoter and exon methylation with gene expression and alternative splicing was assessed using custom R scripts.
The mean FPKM and weighted methylation level were calculated across biological replicates for each sex. The presence of interaction effects in linear models was determined throughout using the anova function in r. Post hoc testing of fixed factors was conducted using the glht function from the multcomp v1.4-12 r package with correction for multiple testing using the single-step method (Hothorn et al., 2008). Correlations were calculated using Spearman's rank correlation rho.

| Additional genome annotation and GO enrichment
Putative promoter regions were defined as 2000 bp upstream from each gene. We chose this size based on previous research (Sinha et al., 2006;Tian et al., 2018). We excluded promoters which overlap any other genomic feature using bedtools (Quinlan & Hall, 2010).
Intergenic regions were determined as regions between the end of one gene and the beginning of the next gene's promoter, excluding any annotated TEs. Gene ontology (GO) enrichment was carried out using the hypergeometric test with Benjamini-Hochberg correction for multiple testing (Benjamini & Hochberg, 1995), using the gostats R package (Falcon & Gentleman, 2007). GO biological process terms were classed as over-represented if they had a q-value <0.05. revigo (Supek et al., 2011) was used to visualize GO terms and obtain GO term descriptions. GO terms for genes with different levels of methylation were tested against a background of all genes. GO terms for genes which show female/male overexpression were tested against a background of all genes identified in the RNA-Seq data. GO terms for genes which show extreme female/male overexpression were tested against a background of all differentially expressed genes.

| Sex-biased gene expression and alternative splicing
All RNA-Seq samples generated between 66.9 and 84.1 million paired-end reads with an average mapping rate of 87% (Appendix S1: 1.0.1). Genes showing different levels and patterns of sex bias are likely subject to different evolutionary processes modulating their expression and sex-specificity (Wang et al., 2015). Therefore, in this study we distinguish three general categories of sex-biased genes. The first category contains sex-biased genes, defined as having >1.5-fold difference in expression between the sexes (q < 0.05).
The second contains extremely sex-biased genes, which are those that show >10-fold difference in expression between the sexes (q < 0.05). The third category consists of sex-limited genes, that is, those with some level of expression in one sex but no detectable expression in the other sex (FPKM =0).
Planococcus citri shows extreme sex-specific expression with many genes showing complete sex-limited expression ( Figure 2a).
We identified a total of 10,548 significant genes with sex-biased expression between P. citri males and females (Figure 2b, Appendix S1: 1.0.2). This is 26.5% of the estimated 39,801 genes in the P. citri genome and 54.7% of all genes identified as expressed in at least one sex in the RNA-Seq data (n = 19,282). Of these sex-biased genes, 10,026 show moderate sex-biased expression (q < 0.05 and >1.5 fold-change) with significantly more showing female-biased expression (5270 compared to 4756, chi-squared goodness of fit: GO term enrichment analysis of sexbiased genes show that both female-and male-biased genes are enriched for core biological processes such as biosynthetic processing and carbohydrate metabolism (Appendix S1: 1.0.3). Additionally, We also identify 168 extremely sex-biased genes (q < 0.05 and >10 fold-change, Appendix S1: 1.0.2), with the majority of these showing extreme male-biased expression (140 compared to 28, chi-squared goodness of fit: χ 2 = 74.667, df =1, p <.001, Figure 2c).
There were only three GO terms enriched for extremely biased male genes; these were 'system process' (GO:0003008), 'sensory perception' (GO:0007600) and 'sensory perception of smell' (GO:0007608). Female P. citri are known to produce pheromones to attract males (Bierl-Leonhardt et al., 1981); therefore, it may be that these extremely male-biased genes are involved in pheromone response. There were no enriched GO terms for genes showing extreme expression bias in females.
Finally, we identify 354 sex-limited genes (q < 0.05 and zero expression in one sex, Appendix S1: 1.0.2) in P. citri. Of these, significantly more are sex-limited to males compared to females (204 compared to 150, chi-squared goodness of fit: χ 2 = 8.2373, df =1, p =.01). GO terms enriched for female sex-limited genes include the following: 'growth' (GO:0040007) and 'anatomical structure development' (GO:0048856) among others (Appendix S1: 1.0.3). GO terms enriched for male sex-limited genes include the following: the same three GO terms mentioned above for extreme sex-biased genes as well as 'proteolysis' (GO:0006508) and some other more general terms (Appendix S1: 1.0.3).
Next, we searched for alternative splicing differences between the sexes. In the current genome annotation, (P. citri v0 mealybug. org), 93.13% of genes are annotated as single isoforms. After filtering out genes which also have low expression in both sexes (<10 FPKM), 1235 genes were tested for alternative splicing. A total of 209 genes were found to be significantly alternatively spliced between the sexes, consisting of 423 isoforms (q < 0.05 and a minimum percentage difference of 25%, Appendices S1: 1.0.4 and S2: section 1.2). The GO terms enriched for alternatively spliced genes are varied, including some related to protein modification (Appendix S1: 1.0.5).
We also manually annotated the DNMT1 gene in P. citri and find it shows higher expression in females compared to males (Appendix S2: section 1.3). To further explore this trend, we checked DNMT1 expression in various life stages and sexes via qPCR (methods can be found in Appendix S2: section 1.3), finding the lowest expression is in 3rd instars, embryos show higher expression and the highest expression is present in adults, although no sex differences were observed, contrary to the results from the RNA-seq data (discussed in Appendix S2: section 1.3.). Finally, it is also worth noting P. citri appears to lack the DNMT3 gene, which is also absent in some other insect species (Bewick et al., 2016).
We next checked to see whether any of the same genes show both sex-biased expression and sex-biased alternative splicing. We found that there was a significant overlap of alternatively spliced genes and genes with sex-specific expression bias (112/209), hypergeometric test, p <.001). The majority of these genes (104/112) also show higher levels of male expression compared to higher female expression (chi-squared goodness of fit: χ 2 = 82.286, df =1, p <.001, Figure 2d). There were no GO terms enriched for female bias and unbiased alternatively spliced genes compared to all alternatively spliced genes as a background. However, male-biased alternatively spliced genes were enriched for metabolic processes and 'proteolysis' (GO:0006508; Appendix S1: 1.0.5).
The sex-determination system in P. citri is unknown and alternative splicing of the doublesex gene has been implicated in sexdetermination in the vast majority of insect species (Wexler et al., 2019). We therefore checked the genome annotation to identify any genes orthologous to the Drosophila melanogaster doublesex gene and whether they were alternatively spliced. There are five genes in the current annotation (P. citri v0 mealybug.org) which are anno- We checked the list of differentially alternative spliced genes between sexes and none of these genes above found in P. citri were differentially alternatively spliced, suggesting the method of sex differentiation in this species is not via alternative splicing of doublesex.
We also checked for sex-biased expression of these genes and only two were expressed, g2969 and g36454, the former shows unbiased expression and the latter shows male-biased expression although overall expression levels are low. Finally, it is also worth noting transformer, a gene required for doublesex splicing (Wexler et al., 2019) is not present in the P. citri genome (Appendix S2: section 1.4).
In a CpG context, females have significantly lower methylation levels compared to males, 7.8% ±0.35% and 9.28% ±0.26%, respectively Males and females also show significant differences in the distribution of CpG methylation levels across genomic features (two-way ANOVA, interaction between genomic feature and sex; F = 316.54, p <.001, Figure 3c). As we have previously shown using the female data in Lewis et al. (2020), P. citri females have unusual patterns of DNA methylation compared to other insect species; we confirm that promoters, exons 1-3 and transposable elements (TEs) have significantly higher levels of DNA methylation compared to exons 4+ and introns (Figure 3c, Appendix S1: 1.0.8). We have found this is also the case for males; however, males show significantly higher levels of methylation than females in all features except for promoters (Appendix S1: 1.0.8). Additionally, we also found high levels of intergenic methylation in both sexes, which, along with promoter and TE methylation, is highly unusual in insect species (Bewick et al., 2019).
In order to determine the distribution of methylation levels across

| Sex-specific DNA methylation: gene level
There were 3,660,906 CpG sites found in all replicates with a minimum coverage of 10×. 75.8% of these were classed as methylated in at least one sample by a binomial test and were then used for differential methylation analysis. A total of 182,985 CpGs were classed as differentially methylated between males and females (q < 0.01 and a minimum percentage difference of 15%), which is around 5% of all CpGs in the genome. The majority of these sites are located in exons 1-3, promoters and TEs (Figure 4a). Exons 1-3 have the highest density of differentially methylated CpGs, followed by promoters and then TEs (Appendix S2: section 1.8.).
Due to the unusual occurrence of promoter methylation in P. citri, we investigated sex-specific differences in promoter methylation and exon 1-3 methylation, separately. We find 2709 genes with a differentially methylated promoter (minimum three differentially methylated CpGs and a minimum overall weighted methylation difference of 15%) between males and females and 2,736 genes with differentially methylated exons 1-3 (Appendix S1: 1.1.0). We chose to further examine regions with a minimum of three differentially methylated CpGs to reduce background noise from possible false-positive calls. A total of 8400 promotor regions contain at least one differentially methylated CpG and 6177 exon 1-3 regions.
A significantly higher number of genes with differential promoter methylation were hypermethylated in females compared to males (2645 in females and 64 in males, chi-squared goodness of fit, χ 2 = 2459, df =1, p <.001, Figure 4b,c). This was also the case for genes with differential exon methylation, with 2709 hypermethylated in females and 33 hypermethylated in males (chi-squared goodness of fit, χ 2 = 2611.6, df =1, p <.001, Figure 4b,d). In females, there is also a significant overlap of genes showing both hypermethylation of the promoter region and exons 1-3 (hypergeometric test, p <.001, Figure 4b).
As males show mostly low-to-medium levels of methylation throughout the genome, we analysed the distribution of methylation levels for each feature determined as differentially methylated.
We found that the average level of methylation in males for male hypermethylated promoters is 0.12 ± 0.07 (mean ± standard deviation) and for exons is 0.14 ± 0.12 (Appendix S2: section 1.8), meaning the minimum 15% threshold difference applied translates into a small actual difference in methylation between males and females. sites should be interpreted with care. We also carried out a GO enrichment analysis for differentially methylated genes but as these genes appear to have no direct effect on gene expression (described later), we discuss the results of this and the relevance of such an analysis in Appendix S2: section 1.9. The same relationship is found between gene expression and methylation of exon 1-3 (Appendix S2: section 1.10).

| Relationship of DNA methylation and expression
On a genome-wide scale, the relationship between gene expression and methylation becomes more apparent (Figure 5c). Genes with no promoter methylation and low levels of promoter methylation have significantly higher expression levels compared to genes with medium and high promoter methylation (linear model: low methylation bin: df =63930, t = 4.93, p <.001, no methylation bin: df =63930, t = 4.047, p <.001, Figure 5d). Again, there is no interaction between sex and methylation bin (two-way ANOVA: F 4,7 = 0.998, p =.392).
The results for exon 1-3 methylation are similar; however, only the low methylation bin has significantly higher expression than genes with medium, high or no exons 1-3 methylation (Appendix S2: section 1.10).

| Relationship of differential DNA methylation and differential expression
If DNA methylation is a causative driver of changes in gene expression, we would expect that differentially methylated genes between sexes are also differentially expressed. Given that higher

(d)
methylation is associated with lower expression in this species, we would also expect that down-regulation of gene expression is associated with higher methylation. However, on a single-gene level, we found there is no clear relationship between the level of differential promoter methylation and the level of differential expression of the corresponding gene (Figure 6a). This is also the case for exon 1-3 methylation (Appendix S2: section 1.11).
Additionally, genes that are hypermethylated in female promoter regions are enriched for genes that show significant expression bias in both females (overlapping genes =113, hypergeometric test with Bonferroni correction, p <.001) and males (overlapping genes =92, hypergeometric test with Bonferroni correction, p =.024). Genes that are hypermethylated in female exons 1-3 are enriched for genes with just female-biased expression, the opposite of our prediction (overlapping genes =138, hypergeometric test with Bonferroni correction, p <.001, Appendix S2: section 1.11). Finally, male hypermethylated genes are not significantly enriched for any genes which show sex-biased expression but male hypermethylated promoters are enriched for unbiased genes (overlapping genes =14, hypergeometric test with Bonferroni correction, p =.021, Appendix S2: F I G U R E 5 (a and b) Scatter graphs of expression levels of every gene plotted against the mean-weighted methylation level across replicates of each gene's promoter region for males and females, respectively. Each point represents one gene. The lines are fitted linear regression with the grey areas indicating 95% confidence intervals. (c) Genes were binned by mean-weighted methylation level of the promoter region across replicates and the mean expression level of each bin has been plotted for males and females. The lines are LOESS regression lines with the grey areas indicating 95% confidence areas. (d) Violin plots showing the distribution of the data via a mirrored density plot, meaning the widest part of the plots represent the most genes. Weighted methylation level per promoter per sex, averaged across replicates, was binned into four categories, no methylation, low (>0-0.3), medium (0.3-0.7) and high (0.7-1). The red dot indicates the mean with 95% confidence intervals We then assessed the overall methylation levels of differentially expressed genes. We found the average promoter methylation level of differentially expressed genes is higher than for unbiased genes in both sexes (linear model: df =27375, t = −10.136, p <.001, Figure 6c).
The same differences are also observed with exon 1-3 methylation (Appendix S2: section 1.11). We then checked to see if a specific set of sex-biased genes, such as those which are sex-limited, drive this overall methylation difference observed. We found no significant difference in methylation between biased, extremely biased and sex-limited categories. We also found the pattern of higher male promoter/exon 1-3 methylation compared to female methylation is the same in most cases (Appendix S2: section 1.11). Finally, it is worth noting we also found annotated genes which were not present in the RNA-Seq data set had considerably higher methylation levels in males and females compared to genes which were expressed in either sex (Appendix S2: section 1.11).

| Relationship of DNA methylation and alternative splicing
Exonic DNA methylation has been associated with alternative splicing in some insect species (Bonasio et al., 2012;Li-Byarlay et al., 2013;Marshall et al., 2019). Therefore, we tested for a relationship between DNA methylation and sex-specific alternative splicing in P. citri. We found that unlike differentially expressed genes, the promoter methylation levels of alternatively spliced genes are lower than nonalternatively spliced genes (linear model: df =27,612, t = −3.772, p <.001). The same pattern is also observed with exon 1-3 methylation (Appendix S2: section 1.11). Additionally, alternatively spliced genes which also show sex-specific expression bias do not significantly differ in their promoter or exon 1-3 methylation levels compared to alternatively spliced genes which show unbiased expression (Appendix S2: section 1.11).
We also then checked to see if alternatively spliced genes were also differentially methylated between sexes. We found only one significant overlap of genes which are both alternatively spliced and differentially methylated (Appendix S2: section 1.11), a single gene was common between alternatively spliced genes which show male expression bias and genes with male promoter hypermethylation (hypergeometric test with Bonferroni correction, p =.034). However, it is likely this overlap is significant due to the small gene lists rather than due to biological significance.

| DISCUSS ION
In this study, we investigated the relationship between sex-specific gene expression and DNA methylation in the mealybug, Planococcus citri, a species with extreme sexual dimorphism and genomic imprinting (PGE). Our major findings include the following: the identification of vastly different genome-wide methylation profiles between the sexes, high levels of intergenic methylation-especially in males, and no relationship between differentially expressed genes and differentially methylated genes, indicating cis-acting DNA methylation does not regulate sex-specific differences in adult gene expression.
We hypothesize that the DNA methylation patterns we observe can be explained by several mechanisms acting simultaneously: (a) the higher and more even distribution of methylation across the male genome could be a cause or consequence of the heterochromatinization of the paternal genome in males, (b) the regulation of a subset of mostly nonsexually dimorphic genes through promoter/ exon methylation in both sexes, (c) the hypermethylation of certain promoters and exons reducing expression in females, possibly to balance expression level between the sexes as a mechanism of ploidy compensation.

| PGE may explain uniform DNA methylation in males
We have identified extreme sex-specific differences in DNA methylation across the genome of P. citri. Most notably, overall higher genome-wide methylation levels in males manifest as low, uniform levels across the genome in comparison to a more targeted bimodal pattern of DNA methylation in females. To our knowledge, this type of sex-specific pattern has not been reported in any other species to date. We have also confirmed promoter methylation in both sexes, which is highly unusual in insects (Lewis et al., 2020). However, we should note UTR regions have not been annotated in P. citri and UTR methylation has been identified in other insect species, for example the aphid M. persicae (Mathers et al., 2019). UTRs are, however, annotated in Strigamia maritima (a species of centipede) which also displays promoter DNA methylation with similar function to what we and Lewis et al. (2020) find in P. citri, leading us to believe UTR methylation does not fully account for the results we have obtained.
We hypothesize these patterns, along with the identification of intergenic DNA methylation, are a result of the unusual reproductive strategy employed by this species, paternal genome elimination. Males with PGE have approximately half of their genome in a heterochromatic state (Bongiorni & Prantera, 2003;Brown & Nur, 1964;de la Filia et al., 2020;Hughes-Schrader, 1948). In mammals and plants, DNA methylation is associated with the formation of heterochromatin (Suzuki & Bird, 2008). Previous research has found DNA methylation differences between the paternal and maternal chromosomes in mealybug species, although studies do not agree upon which chromosome set shows higher levels of DNA methylation (Bongiorni et al., 1999;Buglia et al., 1999;Mohan & Chandra, 2005). It is therefore likely the differences in the pattern of DNA methylation between the sexes may be driven by the condensed paternal chromosomes in males. Future work utilizing reciprocal crosses to identify parent-of-origin DNA methylation at base-pair resolution throughout the genome would further clarify the role of DNA methylation in chromosome imprinting in this species.
While differences in DNA methylation have been associated with the different parental chromosomes, it is the modifications of histones which have been directly linked to the formation of heterochromatin in P. citri (reviewed in Prantera & Bongiorni, 2011).
Most recently Bain (2019) showed that both the H3K9me3-HP1 and H3K27me3-PRC2 heterochromatin pathways are involved in the condensation of the paternal chromosomes in males.
Additionally, non-CpG methylation is also thought to exist in mealybugs in CpA and CpT contexts (Deobagkar et al., 1982) and the genes coding for the necessary enzymatic machinery for these modifications have recently been identified in the mealybug Maconellicoccus hirsutus (Kohli et al., 2020). Although we did not find methylation levels above 0.2% in any non-CpG context (Appendix S1: 1.0.7). These studies suggest PGE is likely mediated by multiple interactions between a variety of epigenetic mechanisms within the genome.

| DNA methylation in females may be involved in ploidy compensation
Another striking pattern we observe is the hypermethylation  (Glastad et al., 2014). The aphid Myzus persicae also shows male hypermethylation on the X chromosome which appears as a single copy in males (Mathers et al., 2019). Although it should be noted female aphids show much higher DNA methylation in the autosomes which are diploid in both sexes. However, it is known in mammals that DNA methylation serves multiple functions in the genome (e.g. Edwards et al., 2017) and this has also been suggested to be the case in insects with the function of DNA methylation potentially changing depending on the genomic context (Glastad et al., 2018). In the examples noted above, higher methylation has been identified in the sex/chromosome which is in the haploid state. DNA methylation in these species is associated with elevated, stable gene expression (Hunt et al., 2013;Mathers et al., 2019), suggesting methylation in these examples may serve to increase expression levels to compensate for single-gene copies. We find a negative relationship between DNA methylation and gene expression in P. citri, suggesting higher methylation in females may serve to decrease expression of certain genes to mirror the haploid expression levels of males. This is further supported by our finding that female hypermethylated genes show overall similar expression levels in both females and males. To test this idea, the expression levels of nonsex-biased genes from each parental chromosome set in both males and females should be assessed. Balanced expression levels would suggest some form of ploidy compensation.
We find no consistent overlap between differentially methylated genes and differentially expressed genes. This suggests that cisacting DNA methylation is not regulating sex-specific gene expression. However, if DNA methylation does indeed play a role in ploidy compensation, we would expect to see no overlap with differentially expressed genes. These findings further support the idea that DNA methylation is involved in chromosome-wide processes, such as paternal chromosome condensation in males and possibly ploidy compensation in females. Indeed, a recent RNAi study, which knocked down DNMT1 in the mealybug Phenacoccus solenopsis, found phenotypic changes in males and females, with females changing colour and losing their waxy coating and males displaying wing abnormalities (Omar et al., 2019). This supports this idea that DNA methylation is involved in the generation of sex differences in mealybugs.
However, another RNAi study in the Hemipteran, Oncopeltus fasciatus, revealed that depletion of DNA methylation did not result in changes in gene or transposable element expression but did lead to aberrant egg production and follicle development (Bewick et al., 2019). Thus, suggesting a functional role for DNA methylation that is independent to specific gene expression. It is also worth noting that previous work in insects has found conflicting evidence for the role of DNA methylation in differential gene expression. Wang et al. (2015) found no correlation between methylation and sex-specific expression in a species of Nasonia, whereas Mathers et al. (2019) found differentially methylated genes between aphid sexes were Finally, while this study focuses on a possible role for DNA methylation in sex-specific gene expression, it is important to consider alternate functions within this context. It has been suggested DNA methylation may be involved in larger genomic processes such as nucleosome positioning (Lewis et al., 2020) and gene duplication maintenance (Dyson & Goodisman, 2020). Some proportion of DNA methylation may also be related to the underlying genotype (Yagound et al., 2020). While exploring these functions is beyond the scope of this study future work, exploring the role of DNA methylation in P. citri should consider these processes in an allele-specific manner. If there are indeed large differences between chromosomes (particularly in males), then this may cloud any results if each chromosome set is not considered independently.

| Sex-specific expression and splicing mirror extreme sexual dimorphism
In addition to our key findings above, we have also identified sexspecific gene expression and alternative splicing. P. citri have no sex chromosomes meaning that males and females share the same genetic complement (Hughes-Schrader, 1948). Thus, the observed sexual dimorphism exhibited must be a consequence of differences in gene expression and splicing between the sexes. Indeed, we found that 54% of genes show sex-biased expression, including a subset of genes that are extremely sex-biased and sex-limited. We found that both male-and female-biased genes are involved in core biological processes. Sex-limited genes are likely important in the phenotypic sex differences observed in P. citri, including sensory-related malelimited genes that may be involved in mate recognition through pheromones (Bierl-Leonhardt et al., 1981). Nasonia males also show extreme sex-biased expression of pheromone genes (Wang et al., 2015). The large number of differentially expressed genes we have identified reflects the extreme sexual dimorphism shown in this species ( Figure 1).
We also identified differentially alternatively spliced genes between the sexes and found a significant number of these show male-biased expression. Genome-wide sex-specific alternative splicing has also been identified in aphids (Grantham & Brisson, 2018) and other insects (e.g. Glastad et al., 2016;Price et al., 2018;Rago et al., 2020). Specifically, Grantham and Brisson (2018) found that differentially expressed and alternatively spliced genes had similar GO term enrichment and they suggest both mechanisms serve to independently generate phenotypic differences between the sexes.
Given the significant overlap of differentially expressed and differentially alternatively spliced genes we have found here, it may be that P. citri utilizes expression regulation and alternative splicing of many of the same pathways to generate phenotypic sex differences.
Additionally, Gibilisco et al. (2016) have shown male and female Drosophila utilize alternative splicing differently-males increase diversity in their gene expression profiles by expressing more genes and females express less genes but use more alternative transcripts.
In P. citri, we found generally more female-biased genes compared to male-biased genes but more male-biased alternatively spliced genes, showing that P. citri sexes also employ different mechanisms to generate sex-specific phenotypes.
We did not find any genes orthologous to the Drosophila doublesex gene. Of the genes identified which may have a similar function to doublesex, none appear to be alternatively spliced. Alternative splicing of doublesex is ubiquitous in holometabolous insects, whereas male-biased expression rather than alternative splicing has been detected in some crustaceans (Kato et al., 2011;Li et al., 2018) and a mite (Pomerantz & Hoy, 2015), indicating male-biased expres-

| Future considerations
It is important to bear in mind that the differences we describe in this study are found in adult whole-body samples and thus do not capture expression and DNA methylation biases between tissues and developmental stages, which are known to vary greatly (Grath & Parsch, 2016;Harrison et al., 2015). Recently, both sex-specific and developmental stage-specific expression has been identified in other mealybug species: Phenacoccus solenopsis (Omar et al., 2019), Planococcus kraunhiae (Muramatsu et al., 2020) and Maconellicoccus hirsutus (Kohli et al., 2020). With Kohli et al. (2020)

| CON CLUS IONS
Overall, this study has shown striking differences in the DNA methylome of male and female Planococcus citri, unlike any previously described sex-specific differences in insects. It is likely these differences are due to the unusual reproductive strategy of this species, paternal genome elimination. Based on our key finding of a lack of direct association between differential DNA methylation and differential gene expression, paired with recent findings by de la Filia et al. (2020) that show males display mostly haploid gene expression, we hypothesize DNA methylation may play a trans-acting role in ploidy compensation in this species, although this is speculation and requires experimental testing. Finally, we have identified a large number of differentially expressed genes between sexes mirroring the extreme sexual dimorphism exhibited in this species and we have found no evidence for sex-specific alternative splicing of genes with possible doublesex function in P. citri. In addition to these key findings this study lays the groundwork for future research exploring the role of DNA methylation in genomic imprinting in insects as well as experimental validation studies to identify the interactions between multiple epigenomic mechanisms which may lead to such extreme sexual dimorphism and paternal genome elimination in this species.

ACK N OWLED G EM ENTS
We thank Peter Sarkies, Andrew Mongue and Kamil Jaron for valuable advice and discussion. This work was supported by the NERC