Introduction

Epigenetic modifications regulate gene expression in response to environmental changes and play important roles in the development of organisms and a variety of human diseases (Jaenisch and Bird, 2003; Lence et al., 2017). In addition to DNA and chromatin modifications, which are well studied, more than one hundred RNA chemical modifications have been identified in cells to date (Lee et al., 2014; Roundtree et al., 2017a). As a marker of post-transcriptional regulation, RNA modifications participate in almost all aspects of RNA metabolism (Roundtree et al., 2017a). N⁶-methyladenosine (m⁶A) is the most prevalent internal modification in many eukaryotic messenger RNAs (mRNAs) and long noncoding RNAs (lncRNAs) (Lee et al., 2014; Roundtree et al., 2017a; Yang et al., 2018), which widely affects RNA alternative splicing (Dominissini et al., 2012), export (Roundtree et al., 2017b), stability (Wang X et al., 2014), and translation (Meyer et al., 2015). RNA m⁶A modifications have been found to be enriched on the transcripts of genes that regulate development and cell fate specification (Dominissini et al., 2012; Meyer et al., 2012; Geula et al., 2015), and some m⁶A sites are regulated in a tissue- or disease-specific manner (Zaccara et al., 2019). In addition, the abnormal expression of m⁶A components is related to the tumorigenesis, proliferation, and metastasis of many types of cancers (Pinello et al., 2018; Ma et al., 2019). Therefore, it is of great significance to study RNA m⁶A methylation for revealing the mechanisms of gene expression regulation and human diseases.

At present, there are few studies on RNA m⁶A modification in Drosophila, possibly due to its relatively low abundance (m⁶A/A < 0.2%) (Haussmann et al., 2016), and the mutation of some m⁶A component genes will affect their viability and fertility (Hongay and Orr-Weaver, 2011; Haussmann et al., 2016; Kan et al., 2017). However, as a model organism with specific genetic and developmental advantages, Drosophila remains an excellent tool for studying the roles of epigenetic modifications in gene regulation, individual development, and disease process. Several components of the Drosophila m⁶A methyltransferase complex (Ime4, dMettl14, and fl(2)d constitute the core complex, with vir and nito acting as cofactors) and a m⁶A reader protein (Ythdc1) have been identified, all of which have homologues in mammals (Lence et al., 2017). Deletion of the major methyltransferase gene Ime4 or the reader Ythdc1 causes locomotion defects, and the splicing of sex-determining factor Sxl is affected (Haussmann et al., 2016; Lence et al., 2016; Kan et al., 2017). However, homozygous mutations in fl(2)d, vir, and nito were lethal, suggesting that these subunits have important functions other than methylation (Penn et al., 2008; Yan et al., 2015; Kan et al., 2017; Lence et al., 2017). RNA m⁶A modification has an obvious sexual dimorphism in Drosophila, and reduced m⁶A levels severely decreased the survival of females. It is thought to be due to the derepression of msl-2 caused by aberrantly spliced Sxl, which forms the Male Specific Lethal (MSL) complex that associates with the X chromosome (Haussmann et al., 2016).

Dosage compensation is a widespread phenomenon in unbalanced genomes (Lucchesi, 2018; Birchler and Veitia, 2021). The deletion or duplication of some chromosomes rather than the whole chromosome set leads to genomic imbalance, that is, aneuploidy (Orr et al., 2015). Aneuploid variation is usually detrimental to organisms (Birchler and Veitia, 2012; Orr et al., 2015; Birchler and Veitia, 2021), and is associated with developmental abnormalities, mental retardation, and various congenital defects (Williams et al., 2008; Huang et al., 2021; Sanchez-Pavon et al., 2021), possibly due to disorders in their gene expression systems (Prestel et al., 2010; Letourneau et al., 2014). Studies across species have pointed out that there is genome-wide trans modulation in aneuploidy, and the genes on the varied chromosomes are compensated to a certain extent, while genes located on the rest of the genome are mainly regulated in the opposite direction to the changes of chromosome numbers, which is known as the inverse dosage effect (Birchler and Veitia, 2012, 2021; Sun et al., 2013c; Hou et al., 2018; Shi et al., 2021). Histone modification (Zhang et al., 2021a), chromatin remodeling (Birchler, 2016), lncRNAs (Zhang et al., 2023) and microRNAs (Shi et al., 2022) have all been shown to play a role in genomic imbalance.

Because of the haploinsufficiency for X-linked genes, heterogametic individuals in organisms with XY sex determination systems could be regarded as analogous to aneuploidy (Disteche, 2016), including humans and Drosophila. Some studies have linked histone H4 lysine 16 acetylation (H4K16ac) and non-coding roX RNAs to the dosage compensation of Drosophila (Park et al., 2010; Conrad et al., 2012). On the other hand, the compensation of X chromosome in human is thought to be regulated by lncRNA X-inactive specific transcript (XIST) (Jordan et al., 2019). Interestingly, in Drosophila, RNA m⁶A modification indirectly affects gene expression through Sxl and msl-2; while in human, m⁶A methylation of the key lncRNA XIST is necessary for the silencing of gene transcription on one of the female X chromosomes (Patil et al., 2016). Moreover, most tumor cells have genomic instability and high levels of aneuploidy (Ben-David and Amon, 2020; Chiarle, 2021), and at the same time, m⁶A component genes are often aberrantly expressed in various cancers (Pinello et al., 2018; Ma et al., 2019).

An increasing number of studies have found that dosage-related effects in aneuploidy may be the integration of multiple modulations rather than through a single mechanism (Prestel et al., 2010; Birchler, 2016). The effects of genomic imbalance are complicated, and the model of dosage compensation and global gene regulation in Drosophila has been extended to autosomal aneuploidies and sex chromosome aneuploid metafemales where MSL complexes are not assembled (Sun et al., 2013b,c; Zhang et al., 2021b). In addition, genomic imbalance and trans regulatory mechanisms in other species such as maize, Arabidopsis, and humans have also been investigated (Hou et al., 2018; Raznahan et al., 2018; Shi et al., 2021; Yang et al., 2021; San Roman et al., 2023). To reveal the role of RNA m⁶A modification in unbalanced genomes, we studied the dynamic changes and regulatory functions of m⁶A methylation under genetic imbalance conditions using autosomal and sex chromosome aneuploid Drosophila maintained in our laboratory. Meanwhile, dosage-sensitive modifiers, differential alternative splicing events, the MSL complex and other factors that may mediate the relationships between RNA m⁶A modification and the dosage-related effects of aneuploidy were also investigated. In summary, we provided a comprehensive picture of RNA m⁶A methylation in unbalanced genomes.

Results

The responses of m⁶A components under genomic imbalance

RNA m⁶A methylation is a reversible epigenetic modification, and its dynamic process is mediated by m⁶A methyltransferases (writers), demethylases (erasers), and m⁶A recognition proteins (readers) (Lee et al., 2014; Yang et al., 2018; Zaccara et al., 2019) (Figure 1A). We first detected the expression of m⁶A components in Drosophila larvae with karyotypes of normal diploids and trisomies that have an additional chromosome arm using RT-qPCR to determine whether RNA m⁶A dynamics are affected by correlative enzymes in unbalanced genomes (Figure 1B; Figure 1—figure supplement 1A,B). It was found that the transcription levels of most m⁶A writers and the major m⁶A reader are down-regulated in aneuploids compared with their respective sex-corresponding controls (Figure 1B). Unlike the other components, the cofactor vir of methyltransferase complex is up-regulated in trisomy 2L females (Haussmann et al., 2016; Lence et al., 2016; Kan et al., 2017). Because RNA m⁶A modification is enriched in the nervous system of Drosophila, the brains of aneuploid third instar larvae were also used to detect the expression of m⁶A components. As expected, a decreased expression trend of m⁶A components similar to that of the whole larvae was observed (Figure 1—figure supplement 1C).

Figure 1.

Previous studies have pointed out that the m⁶A methylomes have temporal and spatial specificity in the development process of organisms, especially for embryonic development and cell differentiation (Meyer and Jaffrey, 2014; Wang Y et al., 2014). Aberrant expression of m⁶A components may cause defects in embryogenesis and even early embryonic lethality (Zhong et al., 2008; Wang Y et al., 2014; Geula et al., 2015; Zhao et al., 2017). Therefore, we designed probes of m⁶A components (Figure 1—figure supplement 1D,E) to examine the mRNA expression and localization patterns during early embryonic development of aneuploid Drosophila using high-resolution TSA-FISH (Lecuyer et al., 2007; Jandura et al., 2017). The results showed that the subembryonic distribution patterns of the five components of m⁶A methyltransferase complex and one m⁶A reading protein are similar in wildtype and aneuploidies (Figure 1C,D; Figure 1—figure supplement 2A-F). At the blastoderm stage (stage 1-5), the probe signals of m⁶A components are widely distributed in the surface cell layer, yolk plasma, yolk cortex, and show a pattern of basal enrichment (Figure 1—figure supplement 2A-F). For the gastrulae at stage 6-11, the transcripts of m⁶A components are mainly located in head, amnioproctodeal invagination, germ band, and as in previous studies (Lence et al., 2016), and show an enrichment in neuroectoderm (Figure 1C; Figure 1—figure supplement 2A-F). At later stages of embryonic development (stages 12-13 and 14-17), the probes are distributed in brain, ventral nerve cord, midgut, and salivary gland (Figure 1C; Figure 1—figure supplement 2A-F).

By observing the subcellular localization of the probe signals of m⁶A components at higher magnification, it can be found that most of the mRNAs of these genes have nuclei-associated localization patterns (Figure 1D; Figure 1—figure supplement 2A’-F’). In the early and late stages of embryonic development, the probe signals of Ime4, dMettl14, fl(2)d, vir, and Ythdc1 form dense small foci near the nucleus, which is a perinuclear distribution (Figure 1D; Figure 1—figure supplement 2A’-F’). For nito, in addition to the perinuclear signals, there is also an obvious signal of intranuclear localization during early embryogenesis, which is manifested as one or two small foci in the blastoderm nucleus (Figure 1D; Figure 1—figure supplement 2E’).

Although there seems to be no difference in the localization of the transcripts of m⁶A components that we detected in aneuploidy and wildtype embryos, the expression levels of these genes, as determined by relative fluorescence intensity, are significantly changed (Figure 1E,F; Figure 1—figure supplement 2A”-F”). Except for the irregular fluctuations at early embryonic stages, the expression levels of most m⁶A components in aneuploids are lower than those in wildtype at more mature stages (Figure 1E,F; Figure 1—figure supplement 2A”-F”), which is similar to the trend we detected in third instar larvae (Figure 1B; Figure 1—figure supplement 1C). The transcripts of m⁶A methyltransferase and reader protein have specific subembryonic and subcellular localization patterns in the development of Drosophila embryos, which may be a mechanism to regulate cellular functions, and ensure appropriate cell growth and differentiation (Lecuyer et al., 2007).

Genome-wide mapping of RNA m⁶A methylation in aneuploid Drosophila

Subsequently, we detected the global levels of RNA m⁶A methylation in wildtype and aneuploid Drosophila larvae to determine whether m⁶A abundance is altered under the condition of genomic imbalance due to modulation by transmethylases. The overall abundance of m6A was represented by the m⁶A/A ratio in total RNA (Figure 2A). It was found that the m⁶A abundance of wildtype males is higher than that of females, but both are lower than 0.2%, which is consistent with the finding that m⁶A levels in Drosophila are relatively low (Haussmann et al., 2016; Lence et al., 2016). Females with triple chromosome 2 left arms (2L) and metafemales with triple X chromosomes have significantly higher m⁶A abundances than diploid females, whereas there is no significant difference between trisomy 2L males and wildtype males (Figure 2A). Therefore, genomic imbalance can affect the m⁶A methylation status to some extent, and this epigenetic modification is different between males and females. However, the m⁶A abundance in aneuploidies did not follow the expression levels of transmethylases.

Figure 2.

To obtain m⁶A mapping across the transcriptome of aneuploid Drosophila, samples with different karyotypes were used for m⁶A methylated RNA immunoprecipitation sequencing (MeRIP-Seq). By identifying consistent peaks in two biological replicates, approximately ten thousand m⁶A peaks were found in wildtype females and males, whereas there were less than half the number of methylation sites in the three kinds of aneuploidies (Figure 2B). When these m⁶A sites were annotated according to the nearest genes, the changes in the number of m⁶A-marked genes in samples of different genotypes were in accordance with that of m⁶A peaks, in which there are about 2,500 m⁶A-marked genes in wildtype and about 1,000 m⁶A-modified genes in trisomies (Figure 2C). Both the number of m⁶A peaks and the number of m⁶A-modified genes are the least in trisomy 2L females (Figure 2B,C). The changes in the number of RNA methylation sites in aneuploids did not coincide with the changes of overall m⁶A abundance, but instead matched the expression of m⁶A components. We speculate that this may be caused by the nonuniformity and heterogeneity of RNA m⁶A modification, including the tissue specificity, the developmental specificity, the different numbers of m⁶A sites in one transcript, the different proportions of methylated transcripts, et cetera (Meyer et al., 2012; Meyer and Jaffrey, 2014; Zaccara et al., 2019). Analysis of m⁶A-modified reads in wildtype and aneuploidy MeRIP samples showed that the expression levels at m⁶A sites in aneuploidies are significantly higher than that in wildtype (Mann Whitney U test p-values < 0.001; Figure 2D), which confirmed our hypothesis.

These results suggest that RNA m⁶A modification exhibits greater heterogeneity in unbalanced genomes.

We next studied the overall characteristics of m⁶A methylation in aneuploid and control Drosophila. m⁶A peaks are localized along all autosomes and sex chromosomes, and the numbers of peaks are correlated with chromosome length (Figure 2—figure supplement 1A). Methylation sites are widely distributed on 5’UTR (30-50%), 3’UTR (10-15%), promoter (15-20%), and internal exons (20-35%), and the aneuploidies appeared to have a higher proportion of 5’UTR peaks and a lower proportion of exon peaks than the wildtype (Figure 2E). The highest values of density distributions of the length of exons containing m⁶A peaks are around 250 bp, and approximately 90% of the exons are longer than the typical 140 bp (Figure 2—figure supplement 1B). De novo motif analysis of m⁶A methylation peaks in each sample showed that the top-ranked motif is consistent with the conserved m⁶A motif DRACH (D=G/A/U, R=G/A, H=U/A/C; Figure 2F). In addition, enrichment analysis of the known motif DRACH showed that it is significantly enriched in all samples, and more than 97% of the m⁶A sites have this sequence (Figure 2—figure supplement 1C).

There are 370 common genes to which the m⁶A peaks were annotated in all genotypes (Figure 2G). We classified the m⁶A-modified genes in each sample and found that they were mostly protein-coding genes, whereas only about 7% were ncRNAs (Figure 2H). As described in previous studies (Meyer et al., 2012; Zaccara et al., 2019), most methylated RNAs have one m⁶A peak (Figure 2—figure supplement 1D-H). However, there are still some genes whose transcripts can be highly methylated and contain more than five m⁶A sites (Figure 2—figure supplement 1D-I). It was found that transcription factors were significantly enriched in 70 genes with more than 5 m⁶A peaks in both wildtype and aneuploidy (Fisher’s exact test p-value = 7.7e-4; Figure 2I; Figure 2—figure supplement 1I). Thus, specific types of genes may be preferentially targeted by m⁶A modification, such as transcriptional regulators (Kan et al., 2017). Analysis of the functions of m⁶A-modified genes in each sample revealed that they were enriched for a large number of common functions, among which the most significant were those related to morphogenesis, development, and growth (Figure 2—figure supplement 2A). The 141 GO terms shared by the five genotypes were summarized and found to include functions about metabolism, regulation, cellular processes, signaling pathways, behavior, and immune response (Figure 2—figure supplement 2B,C). Through pathway enrichment analysis of m⁶A-modified genes, 9 terms were found to be consistently enriched in all samples, including MAPK signaling pathway, Hippo signaling pathway, TGF-β pathway, etc. (Figure 2—figure supplement 2D-F).

Differentially methylated peaks (DMPs) and their associated genes

We then searched for the differential methylation peaks (DMPs) in autosomal and sex chromosome aneuploidies compared with wildtype Drosophila of their corresponding sex. The number of DMPs and the changing trend of DMP methylation status are shown in the Figures (Figure 3A-D). In trisomy 2L females, the number of upregulated DMPs (1,131) is less than that of downregulated DMPs (1,602). In trisomy 2L males, the number of upregulated DMPs (1,521) is slightly higher than that of downregulated DMPs (1,259). In metafemales, the number of upregulated DMPs (1,393) is higher than that of downregulated DMPs (926) (Figure 3A-D). More DMP-associated genes were found in autosomal aneuploidies than in sex chromosome aneuploidy (trisomy 2L female = 2,148, trisomy 2L male = 2,055, metafemale = 1,771; Figure 3E). Also, the total number of peaks after the differential analysis would be different from the number of m⁶A peaks identified in the previous section, because the analysis software merged and recalculated the peaks (Stark and Brown, 2013).

Figure 3.

DMPs are distributed on all chromosomes (Figure 3—figure supplement 1A). Specifically, the MeRIP-Seq signal intensity at the location of DMPs showed that the m⁶A-marked reads in trisomies are denser than those of wildtypes, especially on the triple chromosomes (Figure 3F-H). The density of reads on chromosome 2L of trisomy 2L Drosophila is higher than that on other chromosomes, while the density of reads on chromosome X of metafemales is significantly higher than that of autosomes (Figure 3F-H). Previous studies on aneuploid Drosophila showed that the expression levels of genes on chromosomes with increased copy numbers were mostly similar to those in diploids, and the expression of genes on unvaried chromosomes were widely inverse regulated (Sun et al., 2013c; Birchler, 2016; Zhang et al., 2023). It can be concluded that the changes in the density of reads at DMP sites are not due to a direct gene dosage effect. Similar results were obtained by counting the methylation status of cis and trans transcripts (Figure 3I-K). Overall, the transcripts of genes located on the varied chromosomes in three trisomies have more up-regulated DMPs than down-regulated DMPs, which is more obvious in metafemales (Figure 3I-K). In trisomy 2L females and males, genes localized on other autosomes have a higher proportion of down-regulated DMPs, while the X-linked trans genes show sexual dimorphism and X chromosome-specific response to genomic imbalance. (Figure 3I,J). For metafemales, the transcripts of trans genes have similar numbers of up-regulated and down-regulated DMPs (Figure 3K). Therefore, cis genes in trisomy generally possessed a higher proportion of up-regulated DMP-associated transcripts, whereas the methylation states of trans genes are variable, depending on the identities of varied chromosomes and their genomic locations. The above results give us reason to speculate that the changes of m⁶A methylation may affect dosage compensation of cis-genes and inverse dosage modulation of trans-genes in aneuploidy.

The distribution of DMPs along gene features includes 5’UTR (30-35%), 3’UTR (20-25%), promoter (10-25%), and internal exon (25-35%), among which the proportion of 3’UTR is higher than that of all m⁶A sites (Figure 3—figure supplement 1B). The length distributions of exons with DMPs are similar to that of all m⁶A-modified exons (Figure 3—figure supplement 1C). De novo motif analysis of the sequences where DMPs are located or enrichment analysis of the known motif DRACH both showed that DMP sites are enriched for the conserved motif of m⁶A (Figure 3—figure supplement 1D,E). There are also similarities between the types of DMPs-associated genes and m⁶A-marked genes (Figure 3—figure supplement 1F). The DMP-associated genes of the three aneuploidies have complex overlapping relationships (Figure 3—figure supplement 1G). Most DMP-associated genes have one differentially methylated site, and only a few genes contained more than three DMPs (Figure 3—figure supplement 1H-J). There are 1,042 common DMP-associated genes in all trisomies, of which 71 were transcription factors and were significantly enriched (Fisher’s exact test p-value = 6.9e-5; Figure 4A,B). We speculate that these differentially methylated transcription factors may play an important role in the regulation of gene expression and development of aneuploidy.

Figure 4.

Functional enrichment analysis showed that genes with up-regulated DMPs in three trisomies were enriched for 153 common functions, while the genes with down-regulated DMPs shared only 6 functions (Figure 4—figure supplement 1A,B). Similarly, genes with up-regulated DMPs were enriched for more consistent pathways in all aneuploids (Figure 4—figure supplement 1C,D).

Most of the functions of DMP-associated genes are the same as those of m⁶A-marked genes, and these genes are also enriched for cell fate commitment, dorsal/ventral pattern formation, sex differentiation, post-transcriptional regulation of gene expression, and other additional GO terms. We found that many DMP-associated genes have functions related to gene expression regulation and dosage compensation, including some components of MSL complex (msl-1, msl-3, lncRNA:roX1) and the major sex-determining factor Sxl (Figure 4C). We also noticed that m⁶A component genes fl(2)d, nito, and Ythdc1 themselves are differentially methylated in aneuploidy (Figure 4C), for example, there was a common significantly upregulated m⁶A peak in the 5’UTR of the transcripts of fl(2)d in three aneuploidies (Figure 4D). These data suggest that m⁶A dynamics in aneuploids may be involved in the expression regulation of unbalanced genomes through various genes with regulatory functions.

Relationships between m⁶A and dosage-related modulation of gene expression in aneuploidy

To explore the relationships between RNA m⁶A modification and gene expression regulation under genomic imbalance, we combined the results of MeRIP-Seq and RNA-seq, and performed a comprehensive analysis. The results showed that about 12%-15% of genes differentially expressed in aneuploid Drosophila also belong to DMPs-associated genes (Figure 5A-C). However, except for trisomy 2L males, the DEGs of the other two trisomies did not show enrichment of differential m⁶A methylation (Fisher’s exact test p-values: trisomy 2L male = 0.012, trisomy 2L female and metafemale > 0.05). There are 55 genes that are both differentially expressed and differentially methylated in all aneuploidies (Figure 5D). By analyzing the functions of genes that are simultaneously differentially expressed and differentially methylated, we obtained a number of consistent biological functions, including post-embryonic animal morphogenesis, cell communication, signal transduction, response to stimulation, and so on (Figure 5E), possibly reflecting the important roles of m⁶A-regulated expression in cell signaling and development of organisms.

Figure 5.

By analyzing the crossover of different groups of m⁶A-modified genes and genes with canonical dosage effect, dosage compensation, and inverse dosage effect, we further revealed the relationships between m⁶A-modified genes and dosage-related effects in aneuploidy (Figure 5F-H). The results showed that in trisomy 2L females, cis dosage compensation genes, trans autosomal dosage effect genes, and trans unchanged genes are significantly enriched in m⁶A group with methylation in both trisomy and control; meanwhile, cis dosage effect genes, other autosomal dosage effect genes, and X-linked dosage effect genes are enriched in the group without methylation in both trisomy and control (Figure 5F). On the contrary, in trisomy 2L males, cis dosage effect genes, trans dosage effect genes, and trans unchanged genes are significantly enriched in m⁶A group whose genes are methylated in both trisomy and control; while the group of genes that are not methylated at all is enriched for cis dosage effect genes, cis dosage compensation genes, trans inverse dosage effect genes, and other autosomal dosage effect genes (Figure 5G). Metafemales performed similarly to trisomy 2L female, with significant enrichment of cis dosage compensation genes, autosomal trans dosage effect genes, and trans unchanged genes in m⁶A group where both trisomy and control genes are methylated; and the group without methylation is enriched for cis dosage effect genes and trans dosage effect genes (Figure 5H). These results indicated that m⁶A-modified genes in aneuploid Drosophila females are mainly related to dosage compensation and inverse dosage effect, while genes not modified by m⁶A are mainly related to direct dosage effect. Male aneuploids did not follow this trend.

Furthermore, the distributions of m⁶A modification sites along gene features in genes with classical dosage-related effects were also studied (Figure 5I-K). We found that in wildtype females, trisomy 2L females, and metafemales, a higher proportion of m⁶A sites on genes with canonical dosage compensation and inverse dosage effect are distributed in the 5’UTR than genes with dosage effects (Figure 5I,K). At the same time, this proportion is higher in trisomy than in wildtype. In contrast to females, trisomy 2L males have a higher proportion of 5’UTR RNA m⁶A modification on dosage effect genes (Figure 5J). The above study again demonstrates the sexual dimorphism of RNA m⁶A modification in response to aneuploidy. In addition, it has been suggested that m⁶A residues in the 5’UTR region may have unique regulatory functions (Luo et al., 2014; Meyer and Jaffrey, 2014). Therefore, we hypothesized that the high level of 5’UTR m⁶A modification on genes with dosage compensation and inverse dosage effect might regulate the down-regulation of these genes.

Dozens of dosage-sensitive modifiers have been identified in Drosophila, whose dosage changes can negatively or positively regulate the gene expression across the whole genome like aneuploidy (Birchler et al., 2001; Birchler and Veitia, 2007). Protein-protein interaction (PPI) networks between dosage-sensitive modifiers and differentially expressed genes (DEGs) were constructed to investigate whether RNA m⁶A modification could participate in the gene regulatory networks in unbalanced genomes through these regulators (Figure 5—figure supplement 1A-C). In the three types of aneuploidies, there are more dosage-sensitive modifiers with up-regulated DMPs than with down-regulated DMPs. The regulators ox and Vha55 without DMPs interact with a large number of significantly up-regulated DEGs in trisomy 2L females and metafemales, whereas these two regulators interact with a small number of down-regulated DEGs in trisomy 2L males. Up-regulated DMP-associated regulators wg, Uba1, ap, Atg1, osa, rdx, and sd in trisomy 2L females mainly interacted with significantly down-regulated DEGs (Figure 5—figure supplement 1A); and up-regulated DMP-associated regulators Uba1, osa, and sd in metafemales also interacted with down-regulated DEGs (Figure 5—figure supplement 1C). Up-regulated DMP-associated regulators wg, Kr-h1, and Trl mainly interacted with up-regulated DEGs in trisomy 2L males (Figure 5—figure supplement 1B). In addition, a preponderance of significantly down-regulated DEGs interacted with dosage-sensitive modifiers is observed on trans chromosomes of all Drosophila trisomies (Figure 5—figure supplement 1A-C).

Alterations of m⁶A may be involved in differential alternative splicing in imbalanced genomes

Next, we analyzed the differential alternative splicing events in aneuploid Drosophila. More than 1,000 differential splicing events have been identified in different aneuploids, and about one-third of them are of the type of skipped exon (SE) (Figure 6—figure supplement 1A-C). The biological functions of differentially spliced transcripts are mainly involved in macromolecular fiber organization, locomotion, and growth (Figure 6—figure supplement 1D), and these genes are enriched in heterogeneous pathways in different aneuploidies (Figure 6—figure supplement 1E). Notably, we found that genes with DMPs are significantly enriched for differential alternative splicing events in three kinds of aneuploid Drosophila (Fisher’s exact test p-values < 0.05; Figure 6A-C). Among the genes whose transcripts are differentially spliced, 27%-32% are also differentially methylated (Figure 6A-C). There are 67 genes with both differential alternative splicing and differential m⁶A methylation in all aneuploidies (Figure 6D). The functions of these genes are similar to those of all differentially spliced genes (Figure 6E), but more consistent KEGG pathways are enriched, including endocytosis, mTOR signaling pathway, Hedgehog signaling pathway, and so on (Figure 6F).

Figure 6.

In all trisomy Drosophila, 10 genes are shown to be differentially expressed, with their transcripts also being differentially spliced and differentially methylated. Three of them are illustrated below [su(Hw), Ppn, and CG13124]. su(Hw) is a component of the gypsy chromatin insulator complex, which is a regulatory element that establish independent domains of transcriptional activity (Roseman et al., 1993). Due to the close relationship between su(Hw) and second-site modifiers (Rabinow and Birchler, 1989), and BEAF-32, which is also an insulator DNA-binding protein, has been proposed as a possible inverse dosage regulator (Gurudatta et al., 2012; Zhang et al., 2021b), We speculated that the transcription factor su(Hw) may also be a dosage-sensitive regulator. The data showed that the transcription levels of su(Hw) are up-regulated in all three aneuploidies, and there is a consistent m⁶A modification site with significantly down-regulated methylation (Figure 6G). At the same time, the transcripts of this gene have a common alternative 5’ splice site (A5SS), and its inclusion levels in trisomies are down-regulated, that is, more short transcript isoforms are generated (Figure 6G). Ppn gene encodes an essential extracellular matrix protein that influences cell rearrangements. The expression level of Ppn and its m⁶A methylation in 5’UTR region are both significantly up-regulated in trisomies (Figure 6—figure supplement 1F). In addition, an alternatively spliced exon is significantly less frequently skipped in all aneuploid Drosophila (Figure 6—figure supplement 1G). The third gene, CG13124, which may be involved in involved in regulation of translational initiation, is up-regulated in aneuploids. Two m⁶A sites in its 5’UTR region are methylated at higher levels in all aneuploid Drosophila (Figure 6—figure supplement 1H). Its transcripts also have a common significantly different exon-skipping event (Figure 6—figure supplement 1I).

These results indicate that there are complicated relationships among RNA m⁶A modification, gene expression, and alternative splicing under the condition of genome imbalance. RNA splicing seems to be more closely related to m⁶A methylation than gene transcription. The m⁶A sites located in 5’UTR show remarkable changes in methylation levels in aneuploid Drosophila, and may be involved in the regulation of some differential alternative splicing events, such as exon skipping.

Interactions between m⁶A and Drosophila Male Specific Lethal (MSL) complex

Previous studies have shown that m⁶A components are involved in regulating the alternative splicing of sex-determining gene Sxl in Drosophila, and the deficiency of m⁶A writers or readers will lead to the reduction of female-specific isoforms of Sxl (Haussmann et al., 2016; Lence et al., 2016; Kan et al., 2017). Sxl is also a direct target of RNA m⁶A modification (Kan et al., 2017). We found that Sxl transcripts are both differentially methylated and differentially spliced in three kinds of aneuploid Drosophila (Figure 7A,B). For trisomy 2L females and metafemales, two common m⁶A peaks are significantly up-regulated in the 5’UTR region of Sxl. However, trisomy 2L males have a significantly down-regulated m⁶A peak in the 5’UTR (Figure 7A). Meanwhile, multiple junctions of Sxl transcripts undergo complicated alternative splicing in aneuploid Drosophila, including SE, A5SS, alternative 3’ splice site (A3SS), and mutually exclusive exons (MXE) types (Figure 7A). By checking the distributions of RNA sequencing reads near the male-specific exon (namely the third exon), it can be observed that there are almost no mapped reads on the third exon in wildtype females, trisomy 2L females, and metafemales, while the reads mapped to the second and fourth exons are highly prevalent (Figure 7B). On the contrary, wildtype males and trisomy 2L males have a substantial number of reads on the third exon, accompanied by a smaller number of reads on the adjacent two exons (Figure 7B). Notably, we found that a small number of RNA-seq reads aligned to the male-specific exon appeared in trisomy 2L females, which was identified as an SE-type differential alternative splicing event (FDR = 7.4e-7; Figure 7B). This variation may be related to the abnormal m⁶A methylation levels under the condition of genomic imbalance.

Figure 7.

The abnormal expression of m⁶A components reduces the survival of female Drosophila, which is thought to be probably caused by the expression of downstream MSL complex (Haussmann et al., 2016). To investigate the interplay between the MSL complex and m⁶A modification, we examined the responses of m⁶A regulators in transgenic Drosophila strains with MSL2 mutation or overexpression (Figure 7—figure supplement 1A-E; Figure 7C,D). The results revealed significant changes in the expression profiles of m⁶A regulators in MSL2 transgenic strains, especially the m⁶A reader protein Ythdc1, which increased tens of fold in MSL2 knockdown and overexpressed Drosophila (Figure 7C,D; Figure 7—figure supplement 1E). We also observed that the trends of Ime4 expression in females and males were the same in MSL2-overexpressed Drosophila, whereas there was obvious sexual dimorphism in MSL2-knockdown samples (Figure 7C,D), which may be due to the ectopic assembly of MSL complex in MSL2-overexpressed females (Zhang et al., 2021a). In addition, we also examined the expression levels of MSL2 when Ythdc1 was knocked down, and found that MSL2 was also significantly increased in females (Figure 7—figure supplement 1F). All these results strongly suggest a potential relationship between MSL2 and Ythdc1. Next, we further compared the overall abundance of m⁶A in MSL2 transgenic and wildtype Drosophila (Figure 7E,F). The results showed that m⁶A levels were significantly decreased with MSL2 knockdown and significantly increased with MSL2 overexpression in females (Figure 7E). However, overexpression of MSL2 failed to exert a discernible effect on m⁶A abundance in males, although m⁶A level decreased with MSL2 depletion (Fig.7F). This difference between the sexes may be due to the endogenous MSL2 in males, but not in females.

In the next, we investigated the expression of m⁶A regulators in aneuploid Drosophila overexpressing MSL2, according the results that the MSL complex could be regulated directly or indirectly by m⁶A modification, and the sexual dimorphism of RNA m⁶A modification in response to aneuploidy. It is found that the transcription levels of m⁶A regulators are significantly up-regulated in the brains of most aneuploid larvae that overexpressed MSL2 (Figure 7G; Figure 7—figure supplement 1G). These results are not completely consistent with the quantitative results when MSL2 was overexpressed in diploids, which may be related to the effect of unbalanced genomes. We also used TSA-FISH to detect the expression and distribution of m⁶A components during embryogenesis in aneuploidies overexpressing MSL2 (Figure 7H; Figure 1—figure supplement 2). The subembryonic and subcellular distributions of mRNAs for m⁶A methyltransferases and reading protein did not appear to be affected by ectopic expression of MSL2 in aneuploid Drosophila embryos (Figure 1—figure supplement 2A-F; Figure 1—figure supplement 2A’-F’). But the relative expression of m⁶A components showed diverse dynamics, among which the levels of the most important methyltransferase Ime4 are significantly up-regulated in autosomal trisomy and sex chromosome trisomy with MSL2 overexpression (Figure 7H; Figure 1—figure supplement 2A”-F”). Considering that unbalanced genomes can affect the expression of MSL complex subunits, and the above results show that there is a close relationship between MSL complex and m⁶A modification, we speculate that unbalanced genomes may influence the expression of m⁶A regulators, possibly through the MSL complex.

Relationship of H4K16Ac with m⁶A modification

As an important component of the MSL complex, MOF is a histone acetyltransferase that specifically acetylates histone H4 at lysine 16 (H4K16Ac), thereby affecting chromatin structure and functions, and activating transcription (Kind et al., 2008; Conrad et al., 2012). Previous studies have elucidated the regulatory role of m⁶A in histone modifications, and histone 3 lysine 36 trimethylation (H3K36me3) also plays a role in recruiting methyltransferase complexes to deposit m⁶A markers on RNA (Wang et al., 2018; Huang et al., 2019; Li et al., 2020). To investigate the functions of MOF-mediated H4K16Ac on RNA m⁶A modification, we analyzed RNA-seq data from Drosophila strains overexpressing MOF (Figure 8A; Figure 7—figure supplement 1H-J). We found obvious changes in the expression of m⁶A regulators in strains overexpressing MOF, with the levels of almost all m⁶A regulators being elevated (Figure 8A). It was also observed that after overexpression of MOF, the increasing trends of m⁶A regulators in females and males were not exactly the same, mostly showing more pronounced increases in females, except for Mettl3 (Ime4) and dMettl14. According to the results of RT-qPCR, the expression of Ythdc1 was decreased in MOF knockdown strains and increased in MOF overexpression strains (Figure 8B; Figure 7—figure supplement 1J); meanwhile, the expression levels of MOF were significantly reduced in Ythdc1 knockdown strains (Figure 8C). These results were further verified by polytene chromosome immunofluorescence experiments (Figure 8D,E). In males, the expression level of MOF was significantly decreased in Ythdc1 knockdown Drosophila; while in Ythdc1 knockdown females, the expression level of MOF was not significantly changed compared with wildtype females (Figure 8D,E). These results suggest that MOF-mediated acetylation modification can have a certain effect on RNA m⁶A methylation in a sexually dimorphic manner, which may be due to the endogenous MSL complex and the imbalance of X chromosome dosage in males.

Figure 8.

Next, we investigated whether the m⁶A reader protein Ythdc1 would also have an effect on the levels of H4K16Ac in Drosophila. Quantitative analysis of polytene chromosome immunofluorescence results in Ythdc1 knockdown Drosophila showed that the changes of H4K16Ac levels also showed sexual dimorphism (Figure 8F,G). In males, knockdown of Ythdc1 led to a decrease in H4K16Ac level; whereas in females, knockdown of Ythdc1 did not affect the level of H4K16Ac (Figure 8F,G). These observed changes in H4K16Ac were consistent with the pattern of changes in MOF protein in Ythdc1 knockdown Drosophila strains, suggesting that m⁶A modification may affect H4K16Ac levels through the mediation of MOF. To further study whether unbalanced genomes are involved in the interaction between m⁶A and histone acetylation modification, we also analyzed RNA sequencing data from trisomy 2L Drosophila strains overexpressing MOF. The data showed that overexpression of MOF in trisomy 2L resulted in significant changes in the expression of m⁶A regulators, with a trend different from that observed in diploids (Figure 8A). These results suggest that genomic imbalance might affect the interaction between MOF and m⁶A regulators to some extent, and the potential mechanisms require further investigation.

Discussion

As an emerging epigenetic modification, RNA m⁶A methylation has been found to be involved in almost all aspects of RNA fate and metabolism (Gilbert et al., 2016; Yang et al., 2018). RNA m⁶A modification is also closely related to the development of organisms and a variety of human diseases (Barbieri et al., 2017; Zhao et al., 2017; Pinello et al., 2018; Ma et al., 2019; Liu et al., 2021; Shafik et al., 2021). However, the roles of m⁶A methylation in development and gene expression of aneuploidy have not been studied yet. Aneuploid variation is usually more detrimental than changes of the entire chromosome set due to genomic imbalance (Birchler and Veitia, 2007, 2012). The global changes of gene regulatory networks in unbalanced genomes involves various epigenetic mechanisms, including histone modification, chromatin remodeling, lncRNAs, microRNAs, etc. (Birchler, 2016; Zhang et al., 2021a, 2023; Shi et al., 2022). This study demonstrated that the expression of m⁶A components was altered under genomic imbalance, leading to dynamic changes in the entire methylome. Potential intermediaries by which m⁶A modification could affect trans regulation and achieve dosage compensation in aneuploid Drosophila were also investigated, such as dosage-sensitive modifiers, alternative splicing events, and the MSL complex.

Our experiments show that the expression levels of most m⁶A component genes are significantly down-regulated in aneuploid Drosophila larvae (Figure 1A; Figure 1—figure supplement 1C). Depletion of m⁶A components interferes with the development of animals and plants, especially leading to impaired self-renewal and differentiation of embryonic stem cells, defects in embryonic development, and even early embryonic lethality (Granadino et al., 1990; Horiuchi et al., 2006; Raffel et al., 2007; Luo et al., 2014; Wang Y et al., 2014; Geula et al., 2015). We demonstrated by TSA-FISH that appropriate temporal and spatial specific distributions of the transcripts for m⁶A components are vital during Drosophila embryogenesis (Figure 1C-F; Figure 1—figure supplement 2). The abnormal expression of m⁶A-related genes may affect the development of aneuploid embryos.

In previous studies, the abundance of m⁶A modification is usually positively correlated with the number of m⁶A peaks and m⁶A-marked genes (Luo et al., 2014; Zhu et al., 2023). However, our results obviously did not conform to this rule, with higher m⁶A abundance and fewer MeRIP-Seq peaks in aneuploids (Figure 2A-C). This reflects the complexity and heterogeneity of m⁶A modification. We suspect that in aneuploidy many RNAs may be lost that are methylated at a low level in wildtype, and possess a higher proportion of highly methylated RNAs. Analysis of the expression levels at each m⁶A site confirmed our hypothesis (Figure 2D). Thus, this phenomenon represents an imbalance in m⁶A methylation caused by aneuploidy.

The distribution of m⁶A sites on gene features, m⁶A consensus motifs, biotypes of m⁶A-marked genes, enriched functions and pathways of methylated genes in aneuploidies are similar to those in wildtype. We also found that genes highly methylated in all genotypes are enriched for transcription factors (Figure 2I). It is consistent with previous studies suggesting that transcription regulatory genes may be preferentially targeted by m⁶A (Kan et al., 2017).

We also analyzed the characteristics of DMPs and DMP-associated genes in trisomies. By observing the distributions of MeRIP-Seq reads around DMP sites, it can be found that the densities of m⁶A-marked reads on chromosome 2L are relatively higher in trisomy 2L females and males (Figure 3F,G). Meanwhile, metafemales with triple X have a higher density of m⁶A-marked reads on chromosome X (Figure 3H). Previous studies have found that there is dosage compensation for cis-genes in autosomal and sex chromosome aneuploid Drosophila, and the expression of most genes on varied chromosomes approaches diploid levels (Sun et al., 2013b,c). Therefore, the up-regulation of m⁶A levels is not directly caused by the increased number of chromosomes. A recent study found that the selective enrichment of m⁶A methylation may play a role as a transcript degradation signal in dosage compensation in mammals (Rücklé et al., 2023). Therefore, we speculate that changes in m⁶A levels in aneuploid Drosophila may affect its dosage compensation and inverse dosage effect by regulating the stability of the transcripts. In addition, more up-regulated DMPs are detected in trisomy 2L males and metafemales, while trisomy 2L females have more down-regulated DMPs (Figure 3A-D). Combined with the fact that the survival rate of trisomy 2L female larvae is relatively lower than that of the other two trisomies, it can be speculated that the regulation of DMPs may affect the survival and development of aneuploid Drosophila.

Gene expression in unbalanced genomes is extensively modulated (Sun et al., 2013b,c; Hou et al., 2018; Raznahan et al., 2018; Shi et al., 2021; Yang et al., 2021; Zhang et al., 2021b; San Roman et al., 2023). Previous studies have analyzed the transcriptome data of autosomal and sex chromosome trisomic Drosophila, and found that most of the genes on the triple chromosomes were compensated, and the ratios of their expression levels to wildtype were approximately 1; meanwhile, the expression of genes on other chromosomes was close to two-thirds of that of the wildtype, which was called an inverse dosage effect (Sun et al., 2013b,c; Zhang et al., 2021b). We investigated the relationships between genes with different m⁶A methylation status and genes modulated by classical dosage-related effects (Figure 5F-H). The results showed that for aneuploid females, genes methylated in both trisomy and control are significantly enriched in dosage compensated cis genes and inverse dosage effect trans genes, and genes not methylated at all are enriched in dosage effect genes. However, in aneuploid males, methylated genes are associated with gene dosage effect (Figure 5F-H). We also found that the proportion of 5’UTR m⁶A peaks on dosage compensation and inverse dosage effect genes are increased in aneuploid females (Figure 5I-K). These results provide evidence of sexual dimorphism in the relationships between RNA m⁶A modification and dosage-related effects. m⁶A located in the 5’UTR generally shows higher tissue specificity and richer dynamic changes, and may have unique regulatory functions (Dominissini et al., 2012; Meyer and Jaffrey, 2014; Gilbert et al., 2016). The increased proportion of m⁶A in the 5’UTR of genes with classical dosage-related effects in aneuploidies suggests that m⁶A modification may be involved in the regulation of dosage-dependent genes.

Dosage-related effects of aneuploidy are thought to be caused by dosage-sensitive genes (Shi et al., 2021; Yang et al., 2021). Some dosage-sensitive regulators have been identified in Drosophila, and changes in the dosage of individual regulatory genes can mimic the effects of aneuploidy in the whole genome (Birchler et al., 2001; Xie and Birchler, 2012; Zhang et al., 2021b). Most dosage-sensitive regulators are transcription factors, signal transduction components, and chromatin proteins, which have in common being members of macromolecular complexes or having multicomponent interactions (Birchler et al., 2001; Birchler and Veitia, 2007). We studied the m⁶A modification of dosage-dependent genes and their protein-protein interaction networks in aneuploid Drosophila (Figure 5—figure supplement 1). The results showed that there are complex interactions between dosage-sensitive regulators and differentially expressed genes. Among them, most of the DEGs located on the unvaried chromosomes are down-regulated, indicating that the dosage-sensitive regulators mainly have negative effects on trans target genes. There are more regulatory genes with up-regulated DMPs than down-regulated DMPs in all three aneuploidies, and many of the regulators with up-regulated methylation are connected with interactors that have down-regulated expression.

We comprehensively analyzed the relationships among RNA m⁶A modification, gene expression levels, and alternative splicing. The data showed that differential m⁶A methylation under genomic imbalance appeared to be more closely associated with differential alternative splicing (Figure 5A-C; Figure 6A-C). This phenomenon is reasonable because the mutation or knockout of Ime4, dMettl14, fl(2)d, and Ythdc1 has been found to affect a large number of alternative splicing events in Drosophila, and fl(2)d itself is thought to encode a splicing factor (Penn et al., 2008; Haussmann et al., 2016; Lence et al., 2016). A small set of transcripts are simultaneously differentially m⁶A methylated, differentially expressed, and differentially alternative spliced in aneuploidies, including the transcription factor su(Hw), which may be a dosage-sensitive regulator (Figure 6G-H). Besides, m⁶A modification in the 5’UTR regions may play a special role in some differential alternative splicing events (Figure 6—figure supplement 1F-I).

RNA m⁶A modification in Drosophila has been shown to be involved in the alternative splicing of Sxl (Lence et al., 2017). The deletion of some m⁶A components in females will result in a reduction of female-specific splicing of Sxl and an increase of the inclusion of the third exon, along with phenotypic sexual transformation (Haussmann et al., 2016; Lence et al., 2016; Kan et al., 2017). We found that Sxl transcripts of aneuploid Drosophila are both differentially m⁶A methylated and differentially spliced (Figure 7A-B). Among multiple alternative splicing events, trisomy 2L females have a higher level of the third exon inclusion compared with wildtype females (Figure 7B), which may be related to the differential methylation at 5’UTR of the Sxl transcripts.

Some studies proposed that the deletion of RNA m⁶A methyltransferase Ime4 harms the survival of female Drosophila because of the insufficient inhibition of msl-2 caused by decreased Sxl levels, which in turn leads to up-regulation of X-linked genes (Haussmann et al., 2016). However, other studies have found that the expression of genes on the X chromosome in females with ectopic expression of MSL2 does not increase twofold, that is, the MSL complex does not directly mediate dosage compensation, and the global inverse dosage effect caused by the imbalance of sex chromosomes is the basis for dosage compensation (Birchler, 1981, 2016; Sun et al., 2013a). Therefore, the lethality of the lack of m⁶A in female Drosophila may not be directly caused by up-regulation of X-linked genes through ectopic assembly of MSL complex, and the specific reasons need to be further studied.

The functions of MSL complex in dynamic regulation of global gene expression in aneuploid genomes has been described (Zhang et al., 2021a). Here, we found that MSL2, a structure component of MSL complex, not only affects the expression of m⁶A regulators, but also influences the overall abundance of m⁶A in Drosophila (Figure 7C-F), proving a close relationship between the MSL complex and RNA m⁶A modification. In addition, under the condition of genomic imbalance, the relative expression of m⁶A components was also changed in larvae and embryos of MSL2-overexpressed trisomies (Figure 7G,H). Another component of the MSL complex, MOF, which mediates histone acetylation and transcriptional activation (Kind et al., 2008; Conrad et al., 2012), is also closely associated with the m⁶A reader Ythdc1. Overexpression of MOF increased the expression of Ythdc1 (Figure 8A,B), and in turn, knockdown of Ythdc1 reduced the expression of MOF and the level of H4K16Ac catalyzed by it (Figure 8C-G). These results demonstrate complicated interactions between RNA m⁶A methylation and the MSL complex in unbalanced genomes, which may affect the gene expression, sexual dimorphism, development and survival of aneuploid Drosophila.

Materials and methods

Key resources table

Drosophila stocks and genetic crosses

The Drosophila strains mentioned in this study were all maintained and crossed in our laboratory. The crossing methods have been described previously (Zhang et al., 2021a). Trisomy chromosome 2 left arm (2L) female and male third instar larvae were obtained from the cross of y; C(2L)dp; F(2R) bw females and Canton S males. The metafemale larvae were obtained from the cross of C(1)DX, ywf/winscy females and Canton S males. Aneuploidies overexpressing MSL2 were generated from the crosses of Drosophila with compound chromosomes and MSL2 homozygotes. The MSL2 transgene strain was constructed and validated in a previous study (Sun et al., 2013a). All Drosophila strains were cultured on cornmeal dextrose medium at 25°C. Genes and chromosomal balancers are described in Flybase (https://flybase.org/).

RNA m⁶A methylation quantification

Total RNA was extracted from Drosophila larvae using TRIzol Reagent (Invitrogen) to detect the abundance of m⁶A. The relative quantification of m⁶A methylation was performed using EpiQuik m⁶A RNA Methylation Quantification Kit (Colorimetric) (Epigentek, NY, USA, Cat # P-9005). Specifically, 80 µl of Binding Solution was first added to each well of the plate, and 200 ng of RNA samples, 2 µl of Negative Control, or 2 µl of diluted Positive Control were added to the designed wells. Subsequently, the plate was incubated at 37°C for 90 min. After washing with Wash Buffer, 50 µl of Capture Antibody, 50 µl of Detection Antibody, and 50 µl of Enhancer Solution were added to each well in order, and wells were emptied before adding a new solution each time.

Finally, 100 µl of Developer Solution and Stop Solution were added to each well away from light, and the absorbance was read on a microplate reader at a wavelength of 450 nm.

m⁶A methylated RNA immunoprecipitation sequencing (MeRIP-Seq)

The MeRIP-Seq service was provided by Cloudseq Biotech Inc. (Shanghai, China), and this technology was developed on the basis of published experimental methods (Meyer et al., 2012). In brief, Drosophila larvae of five genotypes, each with two biological replicates, were collected for sequencing. Ribosomal RNAs were removed from total RNA using Ribo-Zero rRNA Removal Kits (Illumina, USA). Immunoprecipitation of m⁶A RNA was performed using GenSeq® m⁶A RNA IP kit (GenSeq, Shanghai, China). The NEBNext® Ultra II Directional RNA Library Prep kit (New England Biolabs, USA) was used for RNA sequencing library construction. High-throughput sequencing was performed using Illumina NovaSeq 6000 sequencers with the paired-end 150 bp protocol.

Analysis of MeRIP-Seq data

The raw sequencing data was first filtered by Trim Galore (version 0.6.10) (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) to remove adapters and low-quality reads. The quality of the data was then assessed by FastQC (version 0.12.1) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Subsequently, clean reads were aligned to the Drosophila reference genome (Drosophila_melanogaster.BDGP6.32.dna.toplevel.fa, downloaded from the Ensembl database) using HISAT2 (version 2.2.1) (Kim et al., 2019). Next, The R package exomePeak2 (version 1.10.0) (Meng et al., 2013) was used for m⁶A peak calling (the screening criteria were log2FC ≥ 1, RPM.IP ≥ 0.5, and score ≥ 5). Motif analysis was performed using HOMER (version 4.11) (Heinz et al., 2010). Differentially methylated peaks were analyzed by the R package DiffBind (version 3.8.4) (Stark and Brown, 2013), which employs the DESeq2 (Love et al., 2014) method. m⁶A peaks with a p-value of 0.1 or less were considered as differentially methylated peaks (DMPs). ChIPseeker (version 1.34.1) (Yu et al., 2015) was used to annotate the peaks, and plotted some of the figures. The profiles and heatmaps of MeRIP-Seq reads around DMPs were generated by deepTools (version 3.5.4) (Ramírez et al., 2014). The signal distribution on the genome was visualized using the Integrative Genomics Viewer (IGV) software (version 2.12.0) (Thorvaldsdóttir et al., 2013).

RNA sequencing data and analysis

The RNA-seq data of aneuploid Drosophila used in this article were generated in a previous study (Zhang et al., 2023) and can be downloaded from the Gene Expression Omnibus (GEO) database (GSE233534). Genome mapping of sequencing data and gene expression quantification were carried out through HISAT2-StringTie pipeline (Pertea et al., 2015; Kim et al., 2019). Differential alternative splicing events were identified by rMATS (version 4.2.0) program (Shen et al., 2014) (with the parameters of -t paired --readLength 150 --cstat 0.0001 --libType fr-firststrand). The method of generating the ratio distributions of gene expression changes was as described previously (Zhang et al., 2022). The plots were generated using ggplot2 (version 3.4.0) (Wickham et al., 2022) in the R program (version 4.2.2). Differential expression analysis was performed by DESeq2 (version 1.38.2) (Love et al., 2014), and the threshold was set to adjusted p-value ≤ 0.05. The functional enrichment analysis was performed by ClusterProfiler (version 4.6.0) (Wu et al., 2021) based on org.Dm.eg.db (version 3.16.0) from Bioconductor (https://www.bioconductor.org/). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway data was obtained from the network (https://www.kegg.jp/). The list of transcription factors of Drosophila was downloaded from AnimalTFDB (Hu et al., 2019). Protein-protein interaction (PPI) relationships were obtained from the STRING database (Szklarczyk et al., 2019).

Real-time quantitative PCR

Total RNA of Drosophila whole larvae, larval brains, or adult heads was extracted using TRIzol Reagent (Invitrogen), and reverse transcription was done with TransScript one-step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech). The sequences of primers designed for RT-qPCR were listed in Figure 1—figure supplement 1A. The usability of these primers was verified by agarose gel electrophoresis (Figure 1—figure supplement 1B). The real-time PCR was performed with TransStart Tip Green qPCR SuperMix (+ Dye II) (TransGen Biotech) using ABI QuantStudio 6 Flex Real-Time PCR System. Relative quantification of gene expression was determined using the 2^-ΔΔCt method.

Fluorescence in situ hybridization (FISH) of Drosophila embryos

Collection of Drosophila embryos and fluorescence in situ hybridization were performed as previously described (Jandura et al., 2017; Zhang et al., 2021a). The probe primers containing flanking T7 promoter elements were listed in Figure 1—figure supplement 1D. After PCR amplification and in vitro transcription, their products were examined by agarose gel electrophoresis (Figure 1—figure supplement 1E). Digoxygenin (DIG)-labelled antisense RNA probes were hybridized to transcripts of interest, and then detected using a succession of anti-DIG antibody conjugated to biotin, streptavidin conjugated to horseradish peroxidase (HRP) and fluorescently conjugated tyramide. This hierarchical process can greatly enhance the probe signals and is referred to as a tyramide signal amplification (TSA) system. All images were acquired with a Zeiss LSM880 laser confocal fluorescence microscope using ZEN software. The same probes for relative fluorescence intensity analysis were photographed using the same parameters. Fluorescence images were processed and analyzed using Fiji (version 1.53c) (Schindelin et al., 2012).

Polytene chromosomes immunostaining

Salivary gland chromosomes immunostaining was performed as previously described (Zhang et al., 2021a). Briefly speaking, the salivary glands from third-instar larvae were first fixed in 3.7% formaldehyde for 1 min, and then dissociated with 50% acetic acid for 5 min. Polytene chromosomes were treated with the following primary antibodies at a dilution of 1:100: anti-SXL (Developmental Studies Hybridoma Bank, M18-s), anti-MSL2 (Santa Cruz, sc-32459), anti-MOF (Santa Cruz, sc-22351), and anti-H4K16Ac (EMD Millipore, 07-329). Finally, fluorescence-conjugated secondary antibodies (Alexa Fluor 488 and Alexa Fluor 594, Jackson Immuno Research) were used for detection.

Acknowledgements

This work was supported by National Natural Science Foundation of China (Grant No. 31871243 and No. 32070566) to L.S. We thank James A.Birchler for revising the manuscript.

Additional information

Funding

Author contributions

Shuai Zhang: Methodology, Investigation, Data Curation, Formal analysis, Writing - Original Draft, Writing - Review & Editing, Visualization. Xinyu Liu, Methodology, Investigation, Validation, Writing - Original Draft, Writing - Review & Editing, Visualization. Ruixue Wang, Methodology, Investigation, Validation, Writing - Original Draft, Writing - Review & Editing, Visualization. Junhan Wang: Methodology, Investigation, Formal analysis. Ludan Zhang: Methodology, Investigation, Validation. Lin Sun: Conceptualization, Writing - Original Draft, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. All authors have read and approved the final manuscript.

Data availability

All raw and processed sequencing data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE253401. The previously published RNA-seq data used in this study are available at the GEO repository under the accession number GSE233534.