Drosophila genetics

All flies were reared at 25°C and kept on standard medium. The mutant Tet alleles are described in [7, 13]; the wild-type allele used in all experiments is w¹¹¹⁸. The stock utilized to examine Robo2 was robo2^x123/CyO [20]. The material used for all whole-genome analysis was either hand dissected third instar larval brains, or, because some experiments necessitated a large input, dissected anterior parts of larvae including the 3 anterior abdominal segments that contain the brain besides other tissues such as imaginal discs, salivary glands, mouth parts and epidermis. Because Tet is highly expressed in the brain and the nerve cell in discs, but not in the other tissues, we call this the Larval Brain Fraction, LBF. Brains and larvae from wt and Tet-GFP third instar larvae were dissected in cold-PBS supplemented with protease inhibitor, snap frozen on dry ice, and stored at -80°C.

Drosophila S2 cells

Drosophila S2 cells were obtained from Thermo Fisher Scientific (Gibco #R69007) and grown in standard Schneider’s medium as described in manufacturing manual.

Immunohistochemistry and imaging

The following antibodies were used for immunolabelling of late stage embryos and chromosomal preparations: mouse anti-Fas2 (Developmental Studies Hybridoma Bank, DSHB), rabbit anti-HRP (Jackson Immunoresearch), mouse anti-Connectin (DSHB), rabbit anti-dsRED (Invitrogen), rabbit and mouse anti-GFP (Invitrogen), mouse anti-H3K4me3 (Invitrogen). Secondary antibodies were purchased from Invitrogen. DNA was labeled with DAPI (Invitrogen). Embryos were collected and fixed via a formaldehyde/MeOH method [13]. Polytene chromosome preparations and staining were performed as in Karachentsev et al. [21]. Images of the ventral nerve cord were obtained using a Leica SP8 using a 40x Objective. Fas2 and HRP labeled embryos were imaged and typically contained 8–10 hemisegments. Hemisegments were examined for midline crossing and in some instances the presence or integrity of the most lateral Fas2+ longitudinal track. Similar imaging and analysis were performed on Connectin/HRP labeled embryos.

LC-MS/MS for 5mC and 5hmC detection and quantification

Mass spectrometry analysis was performed as described previously [22]. Briefly, 3 μL of 10× buffer (500 mM Tris-HCl, 100 mM NaCl, 10 mM MgCl₂, 10 mM ZnSO₄, pH 7.0), 2 μL (180 units) of S1 nuclease, 2 μL (0.001 units) of venom phosphodiesterase I and 1 μL (30 units) of CAIP were added to 1 μg of mRNA from Drosophila wild type and Tet-deficient larval brains, respectively (in 22 μL of H₂O). The mixture (30 μL) was incubated at 37°C for 4 h. The resulting solution was three times extracted with chloroform. The upper aqueous phase was collected and passed through a solid phase extraction cartridge filled with 50 mg of sorbent of graphitized carbon black to remove the salts. The eluate was then dried with nitrogen at 37°C for subsequent chemical labeling and LC-ESI-MS/MS analysis by an AB 3200 QTRAP mass spectrometer (Applied Biosystems, Foster City, CA, USA).

Embryo and larval Tet ChIP-seq

0-12h embryos were collected, processed, and chromatin was prepared according to Yad et al. [23], except lysates were sonicated on a Covaris S2 sonication device (intensity 8, duty cycle 20%, cycle burst 200) for 30 minutes at 4°C to reach fragments ranging from 150–500 bp and then centrifuged at 20,000g at 4°C for 1 minute. Supernatants were collected and centrifuged again for 15 minute to remove debris. Chromatin samples were then snap frozen in dry ice and stored at -80°C until immunoprecipitation in triplicates. All buffers contained cOmplete EDTA-free protease inhibitor cocktail (Roche).

For the larval brain fraction (LBF), 300 frozen larval heads were thawed on ice and 1 ml of NU-1 buffer (5 mM HEPES-KOH pH 7.9, 5 mM MgCl₂, 0.1 mM EDTA pH 8.0, 0.5 mM EGTA pH 8.0, 350 mM sucrose, 1mM DTT). 1% formaldehyde was added to NU-1 buffer before use. Samples were homogenized immediately at room temperature using Dounce with a loose pestle 30 times without foaming for 15 minutes. Samples were filtered first through BD Falcon Cell Strainer 70 μm (Cat No.352350) followed by 50 μm Falcon (Cat No. 340603). Samples were quenched with freshly prepared 125 mM glycine incubated for 5 minutes at room temperature on a shaker and transferred to ice for 5 minutes. Samples were centrifuged at 4000 g at 4°C for 5 minutes. The pellet was washed twice with 1 ml cold PBS and resuspended in 350 μl chilled sonication buffer (50mM HEPES-KOH pH 7.9, 140 mM NaCl, 1mM EDTA pH 8.0, 1% Triton X-100, 0.1% sodium deoxycholate, 1% SDS) and incubated for 20 minutes at 4°C. Lysates were sonicated as described above and chromatin was stored at—80°C until immunoprecipitation.

Chromatin immunoprecipitation

Chromatin samples were thawed on ice and pre-cleared for 15 minutes by rotation in 25 μl of pre-washed binding control magnetic agarose beads (Chromotek). Chromatin was diluted ten-fold in sonication buffer without SDS. 1% of the diluted lysate was recovered and used as input. Diluted chromatin was incubated with 25 μl of pre-washed GFP-Trap MA beads (Chromotek) and rotated at 4°C overnight. Lysates were washed on magnetic stand with 1 ml each low salt RIPA buffer (140 mM NaCl, 1mM EDTA pH 8.0, 1% Triton X-100, 0.1% sodium deoxycholate, 10mM Tris-HCl pH 8.0) (5 times), high salt RIPA buffer (500 mM NaCl, 1mM EDTA pH 8.0, 1% Triton X-100, 0.1% sodium deoxycholate, 10mM Tris-HCl pH 8.0) (2 times), LiCl buffer (250mM LiCl, 1mM EDTA pH 8.0, 0.5% IGEPAL CA-630, 0.5% sodium deoxycholate, 10mM Tris-HCl pH 8.0) (1 time), TE buffer (10mM Tris-HCl pH 8.0, 1mM EDTA pH 8.0) (1 time). All buffers contained cOmplete EDTA-free protease inhibitor cocktail (Roche).

ChIP DNA was eluted by shaking 2 hours at 37°C with 100 μl of elution buffer (1% SDS, 50mM NaHCO₃, 10μg/ml RNaseA), then 4 hours with 0.2μg/ml proteinase K. Beads were concentrated on magnet and elute was recovered. Samples were de-crosslinked overnight at 65°C. Inputs were processed like ChIP samples. DNA was purified by phenol/chloroform/isoamyl alcohol followed by SPRI select beads (Beckman Coulter) and DNA concentration was measured with Qubit fluorometer (Thermo Fisher).

Embryo Tet ChIP-seq library preparation and sequencing

NGS Libraries were made from eluted DNA using the NEBNext Ultra II DNA Library Prep kit (New England Biolabs) according to the manufacturer’s protocol. Briefly, 20 ng of DNA fragments were end-repaired and the blunt, phosphorylated ends were treated with Klenow DNA polymerase and dATP to yield a 3′ A base overhang for ligation of Illumina adapters. After adapter ligation, DNA was PCR amplified with indexed primer for 12 cycles. Libraries were size-selected using Ampure XP beads (Beckman Coulter) to remove adapter dimers. DNA was quantified by fluorometry with the Qubit 2.0 (Thermo Scientific) and DNA integrity was assessed with a Fragment Analyzer (Agilent). The libraries were pooled and sequenced on the NextSeq 500 platform using 75 bp single end sequencing according to manufacturer’s protocol using Reagent v.2.5 at the Waksman Institute Genomics Core. Coverage ranged from 30 million to 60 million tags per ChIP-seq sample.

Larva Tet ChIP-seq library preparation and sequencing

ACCEL-NGS® 1S plus DNA library kit was used to prepare indexed libraries from IP and input DNA. Libraries were pooled respecting equimolarity. Sequencing was performed on Illumina MISeq sequencer in 150 bp paired-end reads.

Embryo Tet ChIP-seq data analysis

Raw reads were trimmed using cutadapt v2.0 [24] to remove adapter and low-quality reads. The processed reads were mapped to the Drosophila melanogaster BDGP6 (dm6) reference genome from Ensembl release 88 using the BWA version 0.7.5-r404 for Chip-seq [25]. For analysis, only unique reads with mapping quality >20 were accepted. Further, redundant reads with identical coordinates were filtered out. Aligned reads were processed by Model-based Analysis of ChIP-seq (MACS2) [26] using Input ChIP DNA as control. For peak calling the MACS2 ‘callpeak’ function was used (-p 1e-2 -g 1.2e+08 -B—nomodel–ext size 147 –SPMR) for each replicate vs. control input. Peaks were selected using the following criteria: p-value <10e-5, fold enrichment over control greater than 10 and a minimal number of reads higher than 50. Bedtools (version v2.24.0) [27] was used to identify overlapping peaks in replicates. A sliding window of 50, 100, 150, 200, 250 and 300 bp around the peak summit (base position of maximum enrichment) was used to determine best range for overlapping peaks. The number of overlapping peaks saturated around window size of 250 bp. Thus, for downstream analysis, windows size of 250 bp was used to identify overlapping peaks in replicates. The Integrated Genomics Viewer (IGV) [28] was used for visualization of ChIP-seq data sets. For visualization in IGV, bigwig peak files were generated using “bdgcmp” function in MACS2 with option “-m logFE -p 0.00001”. Peaks were annotated using the "annotatePeaks.pl" feature of HomerTools [29] with default settings and gtf was obtained from of Ensembl dm6 release 88. De novo motif discovery was carried out on all intersecting peaks of Tet ChIP-seq. DNA sequences (FASTA) were generated from chromosome coordinates produced by peak detection and windowing using the BEDTools. De novo motif analysis was performed using MEME-ChIP [30]. Gene ontology (GO) analysis was done using Database for Annotation, Visualization and Integrated Discovery (DAVID) [31, 32]. Binding profile within gene body was generated using deepTools2 with computeMatrix and plotProfile functions [33].

H3K4me3 ChIP-seq public datasets and analysis

Embryo and larva H3K4me3 ChIP-seq data were obtained from the modENCODE project (GEO: GSE16013) [34]. The analysis was carried out from raw data following the same approach described for Tet ChIP-seq. The overlapping of Tet-ChIP seq peaks and H3K4me3 was computed using BEDTools [27].

Larva Tet ChIP-seq data analysis

Tet-Chip sequencing data were pre-processed using the following steps: the raw sequencing data were first analysed with FastQC. Low-complexity reads were removed with the AfterQC tool [35] with default parameters and Trimmomatic [36] with default parameters was used to remove adapter sequences. The resulting fastq data were again analysed with FastQC to ensure that no further processing was needed. Pre-processed reads were then mapped against the Drosophila reference genome (BDGP6.28) with the bowtie2 algorithm [37] using the ensembl reference transcriptome (version 100). Tet-binding peak regions were identified by applying the MACS2 peak-calling tool [26] to immunoprecipitated (IP) samples, using their input counterpart to estimate background noise (q-value < 0.05). It is worth noting that the “expected genome size” MACS2 parameter was set as the Drosophila genome length excluding ‘N’ bases (i.e., 142 573 024 bp), and summit positions were identified using the MACS2 “-call-summits” option. To avoid identifying extremely large peak regions, the peaks were resized to 100 bp on both sides of the identified summit. Binding profile within gene body was generated using deepTools2 with computeMatrix and plotProfile functions [33].

HydroxyMethylated RNA immunoprecipitation sequencing (hMeRIP-seq)

0-12h embryos were collected, immediately frozen on dry ice, and stored at -80°C until RNA purification. The larval brain fraction (LBF), was dissected, immediately frozen on dry ice, and stored at -80°C until RNA isolation. The RNA immunoprecipitation was performed essentially as described in Dominissini et al. [38]. Briefly, total RNA was isolated using RNeasy Maxi Kit (Qiagen). For each sample 1 mg of total RNA (1 μg/μl) was divided into batches of 45μg and incubated at 94°C in fragmentation buffer (100 mM Tris-HCl pH7.0, 100 mM ZnCl₂) for 40 seconds. Fragmented RNA batches were pooled, and ethanol precipitated at -80°C overnight. RNA samples were washed with 75% ethanol and resuspended in RNase-free water. Fragmentation efficiency was checked on a Bioanalyzer RNA chip (Agilent) to ensure the fragment is 200–300 nt. RNA fragments were denatured by heating at 70°C for 5 minutes, then chilled on ice for 5 minutes. For immunoprecipitation, RNA samples were incubated overnight at 4°C with 12.5 μg of anti-5-hmC antibody (Diagenode rat monoclonal MAb-633HMC) or without antibody as negative control in IP buffer (750 mM NaCl, 50 mM Tris-HCl pH7.4, 0.5% IGEPAL CA-630, RNasin 400 U/ml and RVC 2 mM). 60 μl of equilibrated Dynabeads Protein G (Life Technologies) were added to the samples and incubated at 4°C for 2.5 hours. The magnetic stand beads were washed with 1 ml IP buffer for 5 minutes three times. To elute immunoprecipitated RNA, 1 ml TriPure Reagent (Roche) was added, mixed thoroughly, and centrifuged at room temperature for 5 minutes. Aqueous phase was recovered, and equal amount of chloroform was added, vortexed and aqueous phase was collected after centrifugation and ethanol precipitated at -80°C overnight. RNA was resuspended in nuclease free water and used for library preparation. All buffers contained cOmplete EDTA-free protease inhibitor cocktail (Roche).

hMeRIP-seq library preparation and sequencing

Library preparation was done with the TruSeq ChIP Sample Prep Kit (Illumina) after reverse transcription of pulled-down RNA and synthesis of a second strand (NEB) by Next mRNA second strand synthesis module (NEB)). Briefly, 5 to 10 ng dsDNA was subjected to 5’ and 3’ protruding end repair. Then, non-templated adenines were added to the 3’ ends of the blunted DNA fragments. This last step allows ligation of Illumina multiplex adapters. The DNA fragments were then size selected in order to remove all unligated adapters and to sequence 200-300-bp fragments. 18 cycles of PCR were carried out to amplify the library. DNA was quantified by fluorometry with the Qubit 2.0 and DNA integrity was assessed with a 2100 bioanalyzer (Agilent). 6 pM of DNA library spiked with .5% PhiX viral DNA was clustered on cBot (Illumina) and then sequenced on a HiScanSQ module (Illumina).

hMeRIP-seq data analysis

The processed reads were mapped to the reference genome Drosophila melanogaster BDGP6 (dm6) from Ensembl by using Hisat2 (version 2.1.0) for RNA seq and hMeRIP seq [39]. To analyze gene expression, HTSeq framework, version 0.5.3p9, was used to count the aligned reads in genes [40]. Mode “union” and mapping quality cut-off 20 were used for our analysis. Count-table was normalized so that all samples have the same level of total mapped reads. DEseq2 was used to identify differentially expressed genes [41]. Cufflinks v2.2.1 was applied to calculate the rpkm values [42, 43]. A gene was considered as significantly changed when fold change > = 2 or < = -2 and adjusted p value < 0.05. “SplitNCigarReads” function in GATK (version 3.3–0) (https://gatk.broadinstitute.org/) were used to split reads that contain Ns in their cigar string (e.g., spanning splicing events in hMeRIP-seq data). “rmdup” function of samtools (version 1.3.1) were used to remove a duplicate mapping of reads. Then the same peak calling procedure as ChIP seq data analysis was performed to call peaks and calculate peak score of hMeRIP-seq data. The peaks of hMeRIP-seq were selected using P-value < 10e-5. Peaks of hMeRIP-seq were considered as reduced when the normalized hMeRIP-seq signal (peak score) in control samples was at least 1.4-fold change higher than the signal in Tet depleted samples. The fold change and P-value were calculated using “limma” package in R [44].

Western blot

One hundred third instar larval brains from wild type or Tet^null were dissected and immediately frozen on dry ice. Total protein was isolated from these brains using RIPA buffer and 75 ug of the total protein was loaded to each well. Slit antibody (DSHB, C555.6D, Spyros Artavanis-Tsakonas) was used at 1: 200 dilution and Robo2 antibody [45] was used at 1: 1000 dilution. The western blot signals were detected using IRDye 800CW Infrared Dyes conjugated secondary antibody in LICOR Odyssey CLx imaging system. Signals were quantified using LICOR Image Studio Lite software. See S9 Fig for unprocessed western blot exposure.

Statistical information

Statistical analysis was performed using R or GraphPad Prism 9. Statistics were performed using Student’s t-test or chi-square test unless otherwise specified. Error bars are presented as SEM. P-value < 0.05 is the cut-off for statistical significance.

Tet functions as a 5-methylcytosine hydroxylase and modifies polyA⁺ RNA in S2 cells, embryos, and larval brains

Previously we have shown by dot blot analysis in S2 Drosophila cells and larval brains that the 5hmrC modification was primarily found on polyA⁺ RNA and was strongly reduced in Tet knock-down (KD) cells as well as in larval brains from complete loss-of-function animals (Tet^null) [7]. We have confirmed and quantified these results using ultra-high-performance liquid chromatography tandem mass spectrometry (UHPLC-MS/MS). Measurements of 5mrC and 5hmrC abundance in S2 cells indicate that 5hmrC was strongly enriched in polyA⁺ RNA whereas 5mrC was underrepresented in that fraction as compared to total RNA (Fig 1A and 1B). We then examined changes of 5hmrC and 5mrC in RNA isolated from normal and Tet^null larval brains. We found that 5hmrC was decreased about 5-fold in the mutant brains as compared to control (Fig 1C). Moreover, 5mrC was observed to increase almost 3-fold in the absence of Tet function (Fig 1D). Similar results were found in wildtype (wt) and Tet KD embryos (S1A and S1B Fig). These results confirm and extend our previous antibody-based analyses and indicate that Tet is responsible for much of the conversion of 5mrC to 5hmrC in Drosophila mRNA. As in Tet^null the 5hmrC is not completely abolished an additional enzyme that can modify 5mrC must be present, similar to what is observed in vertebrates [7].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. 5hmrC is found in PolyA⁺ RNA and is controlled by Tet as measured by mass spectrometry.

A. 5hmrC in total and polyA⁺ RNA isolated from S2 cells; B. 5mrC in total and polyA⁺ RNA isolated from S2 cells; C. 5hmrC in total RNA isolated from wild-type and Tet^null larval brain; D. 5mrC in total RNA isolated from wild-type and Tet^null larval brain.

https://doi.org/10.1371/journal.pone.0293894.g001

Tet binds DNA preferentially at the transcription start site of target genes

Members of the Tet protein family are known DNA and RNA binding proteins. Moreover, in vertebrates Tet proteins have been shown to bind DNA at promoter regions to regulate gene expression through active DNA demethylation [19, 46]. We sought to identify the genes that are regulated by Drosophila Tet. We began our experiments by determining if Drosophila Tet also binds DNA and mapping the binding sites. We performed ChIP-seq experiments and mapped Tet-binding peaks genome wide using chromatin isolated from the fly line that expresses the Tet-GFP fusion protein under the endogenous promoter [13]. We used two samples from different stages of development: 3^rd instar larval brain and imaginal discs (larval brain fraction, LBF) and 0-12h embryos. Samples were normalized to input chromatin. As negative control we used chromatin from LBF and 0–12 h embryos that do not carry the GFP gene. In both negative controls that lacked GFP we did not obtain enough ChIP material to allow library preparation and sequencing (see methods). Therefore, no background peaks can be assigned.

Bioinformatic analysis of the LBF ChIP-seq results identified 3413 Tet binding peaks distributed on 2240 genes. An example of Tet binding peak profile is shown in Fig 2A. Tet preferentially occupies promoter regions (Fig 2B) and shows the strongest binding to promoter regions (Fig 2C). Analysis of the Tet bound sites identified a highly conserved CG-rich sequence via MEME-ChIP Motif Analysis (Fig 2D and S2C Fig). This motif is similar to that identified from similar studies of Tet1-bound loci in ESCs [46].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

Genome-wide Tet protein binding sites in Drosophila larval brain fractions, Tet-ChIP-seq analysis: A. Representative gene showing Tet binding peak at the promoter. Arrow indicates promoter orientation; B. Genome wide distribution of Tet occupancy in larva brain fraction. The genomic regions (3’UTR, 5’UTR, exons, intergenic, introns, promoter-TSS transcription start sites, and TTS, transcription termination sites) were defined based on RefSeq gene (dm6) annotations; C. Strength of Tet enrichment on fly genome counted as peak score across the gene body plotted from 3413 peaks; D. Genome wide distribution of Tet binding sites displayed as enriched sequence motif among 3413 peaks identified by de novo motif discovery in this study; E. Binding profile of LBF Tet (red) and H3K4me3 (green) within the gene body ± 5kb; F. 36% of Tet occupied genes on various genomic regions overlapped with the H3K4me3 mark; G. Promoter-associated Tet binding peaks on 40% of genes overlap with H3K4me3 marks.

https://doi.org/10.1371/journal.pone.0293894.g002

The composite model of Tet-binding across the coding region illustrates that Tet occupancy is highest near the promoter and gradually decreases until it undergoes a notable drop at the transcription termination sites (TTS). This closely mirrors the profile observed for H3K4me3, an epigenetic mark associated with actively transcribing regions frequently found at transcription start sites [47] (Fig 2E). While 36% of all Tet peaks co-localize with this chromatin modification (H3K4me3, Fig 2F), 40% of the Tet binding sites at the promoter co-localized with the H3K4me3 mark (Fig 2G).

In embryo samples, we detected 5180 Tet-binding peaks associated with 2578 genes. An example of a Tet binding peak profile is shown in Fig 3A. Tet is enriched throughout the gene body and intronic regions (Fig 3B) however the strength of binding is, like in LBF, strongest at promoters (Fig 3C). The Tet-binding profile across the coding regions is similar to that observed in LBF (Fig 3E). Analysis of the DNA sequences bound by Tet protein in embryos uncovered a highest ranking binding motif that shows significant similarity to the larval Tet consensus sequence (Fig 3D and S2 Fig) and, as with the larval ChIP samples, we observe Tet occupancy to be correlated with H3K4me3 binding sites, at promoters (Fig 3E): 42% of all embryonic Tet peaks co-localized with H3K4me3 chromatin modification marks (Fig 3F) and 51% of the promoter binding sites overlapped with H3K4me3 mark (Fig 3G). In both embryos and LBF, Tet binds to approximately the same number of target genes and 30% of Tet’s targets are identical in both tissues (Fig 3H).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

Genome-wide Tet protein binding sites in Drosophila 0–12 hour embryos, Tet-ChIP-seq analysis: A. Representative gene showing Tet binding peak at the promoter. Arrow indicates promoter orientation; B. Genome wide distribution of embryo Tet ChIP-seq peaks in different genomic regions; C. Strength of Tet enrichment on different genomic regions counted as peak score plotted from 5180 peaks; D. Enriched sequence motif among 5180 embryo Tet ChIP-seq peaks identified by de novo motif discovery in this study; E. Binding profile of embryo Tet (red) and H3K4me3 (green) within the gene body ± 5kb; F. 42% of Tet bound genes in embryo have H3K9me3 modification; G. 51% of genes that show binding of Tet to the promoter that overlap with H3K4me3; H. 27% of Tet bound genes in embryo also have Tet binding peaks in larva brain fraction.

https://doi.org/10.1371/journal.pone.0293894.g003

Our ChIP-seq results indicate that Tet binding sites are distributed throughout the physical map of the genome (S2A and S2B Fig). To confirm these results and show that the Tet-DNA binding domain is sufficient to target Tet to DNA, we constructed transgenic flies carrying a Myc-tagged DNA-binding domain of Tet (CxxC) under the control of the heat shock promoter (hsp70-GAL4::UAS-TetCxxCRFPmyc). We expressed the Tet DNA-binding domain by exposing larvae to heat shock and stained salivary glands with anti-Myc and anti-H3K4me3 antibody. As indicated by Chip-seq, Tet showed many bands distributed on all arms of the chromosomes, but virtually no staining of the chromocenter which contains very few genes. H3K4me3 is also present in a distinct binding pattern on all chromosomes, but in contrast to Tet is abundant in the chromocenter and the nucleolus. These staining results agree with our observation that Tet binds to genes on all chromosomes of Drosophila (S2A Fig).

Our Chip-seq experiments were done in embryos and LBF, two tissues at diverse stages of fly development, but in which Tet protein is highly expressed. In both tissues we identified about 2500 genes that showed significant Tet-binding genome-wide. Tet binding characteristics were similar in both tissues in that the most significant Tet-binding peaks, showing strongest binding, were preferentially located at promoters. About 30% of the Tet DNA-binding sites are identical in embryos and LBF. Thus, it appears that only some of the Tet targets are fixed while others show stage-specific variations throughout development.

Identification of Tet-target mRNAs by hMeRIP-seq in fly tissues

We next determined how many of the genes with Tet-binding peaks also showed 5hmrC modifications of their RNA. To do this we mapped Tet-dependent 5hmrC modifications on RNAs transcriptome-wide in the same tissues we used for our Chip-seq analysis. We first performed hMeRIP-seq on total RNA using basically the same approach we used previously in S2 cells [7]. RNAs isolated from wt 0–12 h embryos and from wt and Tet^null Larval Brain Fractions (LBF) were treated with anti-5hmC antibody or immunoglobulin as negative control, and followed by Next Generation Sequencing (NGS, see methods).

In the embryo, we identified 1815 peaks on 1402 mRNAs. A representative 5hmrC peak profile is shown in Fig 4A. The 5hmrC modification is preferentially associated with coding sequences and a comparison to the expected distribution of peaks shows that the distribution of the modification is not random (Fig 4B). Moreover, as the presence of the 5hmrC modification is not proportional to the abundance of the mRNA the modification appears to function broadly within the transcriptome and is not a regulatory modality restricted either to rare or hyperabundant transcripts (Fig 4C). The 5hmrC-associated sequences identified from these experiments revealed a specific UC-rich motif present within these mRNAs that closely resembles the motif observed in S2 cells and mammalian ESCs (Fig 4D and S3 Fig) [7, 19].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4.

Transcriptome-wide distribution of 5hmrC in Drosophila 0–12 h embryo mRNA, hMeRIP-seq: A. Example of gene showing 5hmrC peak distribution. Arrow indicates promoter orientation; B. Distribution of 5hmrC peaks on embryonic transcripts and comparison of actual and predicted peaks according to the type of structural element within the transcript; C. Distribution of all expressed (gray) or 5hmrC enriched (green) transcripts, showing the number of mRNAs as a function of their expression levels in wt embryo; D. Sequence motif identified in within 1815 5hmrC peaks.

https://doi.org/10.1371/journal.pone.0293894.g004

In mRNA from the wild type LBF, we detected 3711 peaks on 1775 transcripts. A representative profile of 5hmrC enriched peaks in wt and Tet^null is shown in Fig 5A. In wt the peaks were distributed across the gene body (Fig 5B) and 5hmrC marks were found to decorate mRNAs independent of their abundance (Fig 5D). Analysis of the peak sequences indicated the modifications were primarily associated with a UC-rich motif highly related to that identified in embryonic samples (Fig 5F). In mRNA from Tet^null larvae we identified 5,374 peaks in 1710 mRNA. Comparison of mRNAs identified in both the wt and Tet^null samples indicate that the distribution of 5hmrC peaks is similar both in the presence and absence of Tet function. However, in the Tet^null samples, 798 of the transcripts identified or 45% had at least one peak that showed a reduction in the 5hmrC modification relative to wild-type (Fig 6C) and the reduction was most pronounced on intronic and coding region peaks (373 out of 445 peaks or 45% and 429 out of 536 peaks or 46%) compared to the peaks found in the UTRs (5’, 150 out of 625 peaks or 19%, and 3’, 106 out of 577 peaks or 16%; Fig 5C). Thus, within a given mRNA transcript some peaks were affected in Tet^null LBF, while others remained unchanged. These results suggest the preference of Tet to modify specific regions of transcripts. In the absence of Tet we find an increased number of peaks on mRNAs, suggesting that Tet may have an inhibitory effect on any additional 5mrC modifying enzyme.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

Transcriptome-wide distribution of 5hmrC in LBF mRNA, hMeRIP-seq: A. Example of gene showing 5hmrC peak distribution. Arrow indicates promoter orientation; B. Distribution of 5hmrC peaks on wt LBF transcripts and comparison of actual and predicted peaks according to the type of structural element within the transcript; C. Distribution of 5hmrC peaks with reduced peak score in Tet^null compared to the peaks found in the wt LBF; note that peaks in the protein coding sequences and introns are significantly more reduced in Tet^null than are the peaks in the 5’ and 3’ UTR. D. Distribution of all expressed (gray) or 5hmrC enriched (green) transcripts, showing the number of mRNAs as a function of their expression levels in wt LBF; E. Distribution of all expressed (gray) or hmrC enriched (green) transcripts, showing the number of mRNAs as a function of their expression levels in Tet^null LBF; F. Sequence motif identified within 3711 5hmrC peaks.

https://doi.org/10.1371/journal.pone.0293894.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. The 5hmrC modified mRNAs.

A. Volcano plot of mRNAs that are increased (green) or decreased (red) relative to wildtype levels in Tet^null LBF preparations; B. Proportion of modified mRNAs in all 9000 wild type transcripts, and in the decreased and increased portions of mRNAs from Tet^null LBF; note the low level of modified transcripts in these two groups of mRNAs; C. The percent of transcripts that show a reduction of 5hmrC modification in Tet^null compared to wt and the percent of transcript that showed both 5hmrC reduction and Tet binding to the corresponding gene; D. GO term analysis of transcripts that show reduction in 5hmrC modification; E. IGV tracks of a representative gene showing the distribution of indicated peaks along the gene body; F. Pathway analysis of neuronal genes shown in D. ChIP-seq, hMeRIP-seq and RNA-seq data are shown in reads per million with the y-axis. Genomic regions with statistically significant enrichment were measured by -log10 (peak P values); P<10⁻⁸) are indicated. The effects of Tet depletion on 5hmC levels are also represented. The Y axis scale is indicated above each track. Blue arrows show reduction in 5hmrC peaks.

https://doi.org/10.1371/journal.pone.0293894.g006

In addition, 37% of the modified mRNA in embryos were also identified in the LBF, while 30% of the larval modified mRNAs were also present in the embryonic fraction (S4C Fig). Taken together these results suggest that Tet targets a distinct cohort of mRNAs in embryos and larval brains and controls specific 5hmrC modifications along transcripts.

RNA levels in wild type and Tet^null larval brains

Our results indicate that Tet binds to the promoter of a subset of possibly actively transcribed genes and controls the 5hmrC modification of their mRNAs. The modification may have an effect on the stability, processing, and/or translation of the transcripts. To determine if there is a link between 5hmrC modification and mature mRNA levels, we sequenced (NGS) RNA isolated from wildtype and Tet^null LBF. We found that out of 9000 total transcripts the levels of 445 were significantly increased and 115 were decreased in Tet^null LBF (Fig 6A). When we compared these mRNAs with the 5hmrC-modified mRNAs present in LBF, we found that 1716 or ~20% of the total transcripts were modified. However, of these modified mRNAs only 15, or 3% were upregulated in Tet^null, and 13 or 11% were decreased (Fig 6A, 6B). This result indicates that the levels of the vast majority of 5hmrC modified mRNAs do not change levels in Tet^null LBF. Thus, the 5hmrC modification of the mRNAs does not appear to generally control the steady state level of transcripts.

Cellular function of genes controlled by Tet

Tet protein is detected in embryos from blastoderm stage onwards and is most strongly expressed in neuronal tissues and also in cardiac and muscle precursor cells. In third instar larvae, the gene is strongly expressed in the brain and neuronal cells in imaginal discs [13]. It was therefore important to assess if our molecular analyses would agree with this expression pattern and if target genes are associated with neuronal functions. We performed Gene Ontology (GO) analyses of the genes identified via ChIP-seq as well as of the genes encoding the 5hmrC-modified mRNAs that were identified in our hMeRIP-seq analyses in the embryo and the LBF (S5A–S5D Fig). The genes identified in both embryonic and larval samples through both ChIP-seq and hMeRIP-seq all show enrichment for genes involved in axon guidance. When we looked at the GO terms of transcripts that showed a reduction of the 5hmrC modification in Tet^null samples, axon guidance genes were highly represented, in fact, GO terms of transcripts showing reduction of the modification in Tet^null samples identified mostly genes associated with neuronal functions (see highlighted genes in Fig 6C). Pathway analysis showed that ROBO receptor signaling is the most enriched pathway in the group of neuronal genes (Fig 6F).

It is striking that in our two very different experimental approaches, ChIP-seq and hMeRIP-seq we identified genes with overlapping functions (S5A–S5D Fig). The importance of our results is also underlined by the observation that of the transcripts that reduction of 5hmrC levels in Tet^null samples, 40% were derived from genes that also have at least one Tet DNA-binding site (Fig 6C). In LBF samples, 43% of all the transcripts that show 5hmrC modification are derived from genes that have been shown to bind Tet (S4A Fig). In embryo samples, 29% of all the transcripts that showed 5hmrC modification are derived from genes that bind Tet (S4B Fig). Further, 29% of modified transcripts in embryos and 37% of modified transcripts in LBF show 5hmrC marks at both developmental stages (S4C Fig). An example of the experimental IGV tracks of all our results for a gene in the larval CNS and the embryo are shown in Fig 6D and S6A Fig, respectively.

These analyses show that Tet-dependent 5hmrC is often found on mRNAs derived from genes that show Tet binding. Notably, close to 50% of transcripts that show a reduction in the 5hmrC mark in Tet^null tissues are derived from Tet-target gene identified by ChIP-seq. However, the levels of these mRNAs are generally unaffected by the loss of Tet suggesting that the 5hmrC modification does not affect steady state level of mRNAs but other aspects of mRNA function such as translation or localization.

Tet target genes

We used the results above to identify Tet-target genes and sought to determine whether the phenotypic effects of the loss of Tet’s activity were derived from its inability to regulate target mRNAs [7, 13]. We looked for genes that are 1. active in the nervous system where Tet is enriched and 2. showed Tet protein binding to DNA, and 3. whose mRNA showed a reduction in 5hmrC in Tet^null animals. Axon guidance genes as a group frequently showed Tet-DNA-binding and 5hmrC mRNA modification by Tet (Fig 6D) and Robo receptor signaling is the most enriched pathway. Among the genes that fulfilled the three criteria were two well-studied genes that function in Robo receptor signaling, robo2 and slit (S7 Fig). The Slit/Robo signaling pathway is required for axonal pathfinding and the bilateral organization of the CNS in both vertebrates and invertebrates [48]. Robo proteins are transmembrane receptors on axonal growth cones for the secreted Slit ligands. Glial cells present at the midline secrete Slit and signaling between Robo and Slit is essential to inhibit midline crossing of axons through commissures via repulsion [49]. Importantly, Slit has previously been implicated as a target of Tet activity in midline glia [16]. We examined axonal pathfinding in the embryonic ventral nerve cord (VNC) and reasoned that if Tet impinges upon the levels of Robo2 and/or Slit, we should observe midline defects in Tet^null animals like those seen in robo2 or slit mutant embryos. Gross CNS commissural structure is maintained in Tet^null embryos (Fig 7B’, HRP), however, examination of neuronal subpopulations within the longitudinal neuropils indicates frequent pathfinding defects. A well described neuronal subpopulation, Fas2+ neurons, exhibit extensive midline crossing of growth cones in these Tet^null embryos (Fig 7B, arrows; S1 Table). Additionally, the most lateral of the Fas2⁺ longitudinal tracks are often incomplete or absent (Fig 7B, 46%-arrowheads). A second subpopulation of neurons expressing Connectin also appears to be altered in Tet^null VNCs and fails to populate one of the longitudinal tracks compared to wild type (S7B Fig; arrows). These phenotypes are strikingly similar to the axonal pathfinding defects seen in robo2 embryos with Tet’s effects being slightly more severe (Fig 7B and 7C and S1 Table) [49]. We sought to determine whether the reduction of Tet-mediated 5hmrC deposition on the robo2 or slit mRNAs resulted in mRNA species with reduced activity or potential for expression. Thus, we examined genetic interactions between Tet and the Slit/Robo signaling pathway in Tet^null embryos lacking one copy of robo2 or slit. We additionally examined Robo1, a gene that is also involved in midline repulsion but is not 5hmrC modified. Decreasing the dose of Robo2 or Robo1 in a Tet^null background has little effect on Fas2+ axonal pathfinding in comparison to Tet^null alone (S1 Table). The failure to see an effect with Robo2 may stem from the observation that the levels of midline crossing in Tet^null embryos exceeds that seen for robo2^null embryos (S1 Table and [50]). However, reducing the gene dose of Slit by half enhances the midline crossing of Fas2+ neurons in Tet^null embryos (S1 Table; Fig 7B and 7E; 48% vs 32% Tet^null), whereas heterozygous slit embryos show midline crossing in < 1% of segments (Fig 7D). Moreover, Tet^null mutant animals appear to be sensitized towards midline crossing in general when lacking full slit function. Notably, the commissures (red arrowheads, Fig 7E’) are poorly defined, likely due to too many axons inappropriately transiting the midline.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Tet regulates the expression of members of the Slit/Robo signaling pathway.

Stage 16/17 embryonic ventral nerve cords immunolabelling a subpopulation of CNS neurons with Fas2 (A-C) and the general neuronal cell surface marker, HRP (A’-C’). A, A’: wild-type; B, B’: Tet^null/Tet^null; C, C’: robo2^x123/robo2^x123; D, D’: sli²/+; E, E’: sli²/+; Tet^null/Tet^null. Examples of midline crossing are indicated by white arrows and malformed lateral Fas2 tracks are noted with white arrowheads. Red arrows in E’ highlight commissural malformations present in Tet^null/Tet^null embryos with reduced slit dosage. Percentage midline crossing is displayed in the overlay panels. F. Western blot showing Slit and Robo2 proteins in wt and Tet^null/Tet^null 3^rd Instar larval brain extracts. GAPDH is the loading control, S9 Fig for full western; F’. Normalized levels of Slit and Robo2 quantitated via optical densitometry.

https://doi.org/10.1371/journal.pone.0293894.g007

Given that robo2 or slit encode mRNAs that carry 5hmrC mark and exhibit a reduction in a Tet^null background while maintaining normal steady state mRNA levels, we expected Tet to potentially control their protein levels (S7 Fig). Indeed, both proteins were clearly reduced in brain extracts from Tet^null larvae relative to wt (Fig 7F and 7F’). These results support the idea that one function of Tet-dependent 5hmrC modification is to control high levels of translation of specific target mRNAs and that in the context of embryonic axonal pathfinding Tet provides an additional, novel layer of regulation of the medically important Slit/Robo pathway.

Based on all our results we suggest the model shown in Fig 8, we propose that Tet binds, possibly as a complex to DNA binding sites mediated by its DNA-binding domain. The Tet binding sites are preferentially located at promoter regions of genes that also show H3K4me3, generally accepted as a mark of active transcription. We further postulate that Tet binds nascent mRNA through its RNA binding domain or possibly in cooperation with associated proteins (RNA-binding proteins, and with a so far unidentified RNA methyltransferase) to set the 5hmC mark. The 5hmrC marked mRNAs are then exported from the nucleus and recognized by a reader protein that will control the efficient loading of the modified mRNAs onto polysomes, where the mRNAs are proficiently translated. It remains to be determined if the 5hmrC modification is the first step in a demethylation cascade, as 5hmC is in vertebrate DNA.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. The proposed model of Tet functions in RNA modification (see text for description).

https://doi.org/10.1371/journal.pone.0293894.g008

While several aspects of this model need to be investigated our results provide a consistent framework of how Tet and Tet-dependent RNA modifications may function in controlling gene expression. Recently, mutations in human Tet3 have been shown to cause neurodevelopmental delays. It will be interesting to investigated if 5hmrC RNA modification is deficient in the affected patients [51].

Tet controls the 5hmrC modification on mRNA

In mass spectrometry experiments we determined that 5hmrC is strongly enriched in polyA⁺ RNA confirming our previous dot blot results. Because Tet is expressed in Drosophila almost exclusively in nerve cells, we determined the levels of 5mrC and 5hmrC in two tissues that show high Tet expression, wild type 0–12 h embryos and in larval brains. We found that 5hmrC levels are similar to the 5mrC levels in wt embryos and larvae and that in both tissues the level of 5hmrC is decreased about 5 times, while the level of 5mrC is increased by about 3 times. In the absence of Tet and therefore most of the mRNA 5hmC modification an increase in its substrate 5mrC is expected. Together all of our mass spectrometry results are consistent with Tet controlling the transition from 5mrC to 5hmrC, mostly in polyA+ RNA.

In all of our experiments we find that even in the absence of Tet ~20% of the wild type 5hmrC levels remain. This observation suggests the presence of additional enzymes that can modify 5mrC in Drosophila. The existence of additional enzyme(s) contributing to mRNA hydroxymethylation has also been postulated in mouse ESCs [19].

Our mass spectrometry findings and the results from our hMeRIP-seq experiments on larval brain fractions (LBF) and embryos are consistent with what has been previously reported for Drosophila tissue culture cells and for ESCs (Figs 1, 4, 5 and S1, S3 Figs) [19]. We identified ~3000 5hmrC peaks in ~1500 transcripts in S2 cells [7]. In ESCs the number of peaks was 1633 in 795 transcripts [19]. In our in vivo experiments we identified 1815 peaks in 1402 transcripts in embryos, and 3711 peaks on 1776 transcripts in LBF. Of the modified transcripts in embryos 37% were also identified as modified transcripts in the LBF. In all samples the modification peaks centered around a UC-rich consensus motif (S3 Fig). The consistency of the mapping results of the 5hmrC modifications in Drosophila tissue culture cells, embryos, larval brain fraction, and ESCs underlines the probable conserved function of Tet across the species.

The 5hmrC peaks on mRNAs derived from LBF are distributed all along the transcripts, the UTRs, the coding region, and introns. However, in Tet^null LBF peaks in the CDS and introns are significantly more strongly reduced than peaks in the UTRs (Fig 5C). This observation suggests that Tet may target coding sequences and introns specifically. We do not yet understand if modifications in different parts of the transcripts have diverse functions and if they may be controlled by additional enzyme(s).

Drosophila Tet’s DNA binding activity

We found that in both embryos and in LBFs, Tet recognizes a DNA motif similar to the motif bound by Tet1 in vertebrate ESCs (S2C Fig) [46, 53]. A majority of these peaks are associated with coding regions and are frequently found at the promoter. Almost 50% of the peaks overlap with the H3K4me3 mark, an indication that the genes are actively transcribed. The distribution of Tet-binding peaks and the overlap with the H3K4me3 mark agree well with the localization of the Tet-DNA-binding domain on salivary gland chromosomes confirming that the binding sites are found almost exclusively in euchromatin and are distributed on all 4 chromosomes (S2A Fig).

We propose that the selection of target RNAs modified by Tet is at least in part facilitated by Tet’s DNA-binding of specific genes. The concurrence of Tet-DNA binding peaks on genes that also showed Tet-dependent 5hmrC modifications of their mRNA is consistent with this idea. The majority of the genes that show Tet binding and modified mRNAs are divergent in both tissues indicating that in addition to a conserved function of Tet in different neuronal cells, Tet also has a tissue-specific or possibly even cell-type-specific function.

Identifying Tet target mRNAs

Tet is highly expressed in nervous tissues and the loss of Tet function leads to abnormal neuronal functions such as defects in larval locomotion or abnormalities in the circadian rhythm. [13] Our immunoprecipitation of 5hmrC-modified RNAs identified 1775 genes in larval brain fractions. 45% (798) showed a significant decrease in the overall 5hmrC peaks in a Tet^null background. Of the genes with reduced 5hmrC marks, 44% showed Tet-DNA binding. Notably, the mRNAs in which the reduction of the 5hmrC mark was seen were mostly associated with genes that function in different aspects of nerve cell development. First among them are axon outgrowth genes that were also identified in the GO-term analysis as abundant gene categories associated with Tet binding sites and mRNAs carrying the 5hmrC mark (Fig 6D and S5 Fig).

Our initial examination of the developing embryonic ventral nerve cord (VNC) in Tet mutants identified subtle defects in CNS patterning. We then examined subsets of VNC neurons using antibodies to Fas2 and Connectin (Fig 7B, 7B’ and S7B, S7B’ Fig) guided by our molecular results. Overall commissural structure is maintained in Tet^null embryos, however neurons expressing Fas2 show a failure of the midline to repel axon crossing effectively. And so, we looked among the Tet mRNA targets with known functions in axon guidance and found that both slit and robo2 mRNAs were represented. Both genes have Tet-binding sites near the TSS, their mRNA is modified, and the modification is reduced in Tet^null LBF, while their mRNA levels are not significantly changed (S7 Fig). Comparison of the CNS in Tet^null and robo2^null embryos identified a set of overlapping phenotypes with high frequency midline crossing defects of Fas2+ neurons, as well as discontinuities in the most lateral, longitudinal Fas2 and Connectin axonal tracts (for description of embryonic nerve cord see [49]). Notably, these tracts correspond to neurons which express the Robo2 protein [54, 55].

The overlapping phenotypes of Tet, robo2 and slit, together with the molecular data that identified Robo2 and Slit as Tet targets, prompted us to investigate if Robo2 and Slit protein expression was affected by the loss of Tet. Indeed, in Western blots from Tet^null larval brain extracts both Robo2 and Slit protein levels were strongly reduced (Fig 7F and 7F’), indicating that Tet’s profound consequences on VNC patterning occurs, at least in part through the control of expression of the Robo2 and Slit proteins. As Robo2 and slit mRNA levels are not changed in Tet^null LBF (S8 Fig), we suggest that the Tet-dependent 5hmrC modification positively controls the level of translation of the two mRNAs. While we have not investigated the protein levels of additional Tet-targets, we expect that Tet controls protein levels through the 5hmC modification of many target mRNAs. Which step in RNA processing leading to mRNA translation is affected in Tet^null animals will have to be elaborated. Based on our previous results, that showed that 5hmrC modified RNAs are found on polysomes, at least one possibility is that the 5hmrC modification facilitates the loading of the mRNAs on ribosomes [7].

Our work supports a function of Tet in controlling the 5hmrC modification of specific neuronal mRNAs, essential for maintaining translation levels necessary for normal neuronal function, thus adding an additional level of control of gene expression. However, we cannot exclude that Tet has additional functions in controlling gene expression in Drosophila.