The role of Dmnt1 during spermatogenesis of the insect Oncopeltus fasciatus
Epigenetics & Chromatin volume 16, Article number: 28 (2023)
The function of DNA methyltransferase genes of insects is a puzzle, because an association between gene expression and methylation is not universal for insects. If the genes normally involved in cytosine methylation are not influencing gene expression, what might be their role? We previously demonstrated that gametogenesis of Oncopeltus fasciatus is interrupted at meiosis following knockdown of DNA methyltransferase 1 (Dnmt1) and this is unrelated to changes in levels of cytosine methylation. Here, using transcriptomics, we tested the hypothesis that Dmnt1 is a part of the meiotic gene pathway. Testes, which almost exclusively contain gametes at varying stages of development, were sampled at 7 days and 14 days following knockdown of Dmnt1 using RNAi.
Using microscopy, we found actively dividing spermatocysts were reduced at both timepoints. However, as with other studies, we saw Dnmt1 knockdown resulted in condensed nuclei after mitosis–meiosis transition, and then cellular arrest. We found limited support for a functional role for Dnmt1 in our predicted cell cycle and meiotic pathways. An examination of a priori Gene Ontology terms showed no enrichment for meiosis. We then used the full data set to reveal further candidate pathways influenced by Dnmt1 for further hypotheses. Very few genes were differentially expressed at 7 days, but nearly half of all transcribed genes were differentially expressed at 14 days. We found no strong candidate pathways for how Dnmt1 knockdown was achieving its effect through Gene Ontology term overrepresentation analysis.
We, therefore, suggest that Dmnt1 plays a role in chromosome dynamics based on our observations of condensed nuclei and cellular arrest with no specific molecular pathways disrupted.
Cytosine methylation of insect DNA remains an enigma. The addition of a methyl group to the nucleotide cytosine is important to cell function and survival of vertebrates and plants [27, 42]. For many domains of life, cytosine methylation is an epigenetic mechanism that negatively regulates gene expression . The cytosine methylation of insects is different. Cytosine methylation is highest within exons and not gene promoter regions [14, 17] Gene body methylation levels across insects is associated with highly and broadly expressed genes [13, 17, 28, 49]. However, differences of cytosine methylation are rarely directly associated with difference of gene expression [13, 34], but might be associated with reduced variation of expression . The function of methylation for insects appears to be lineage specific [32, 47] and direct functional manipulation of Dmnt1 implicates a role in reproduction independent of cytosine methylation [4, 8, 22, 48]. In addition, cytosine methylation itself and its machinery are variably conserved across most insect lineages [7, 13, 14, 39], although most insects have copies of methyltransferases . Indeed, this variation can be extreme. Some insects, such as Drosophila melanogaster, lack DNA methyltransferases and others, such as Tribolium castaneum, have Dmnt1 but lack cytosine methylation in their genome [7, 39]. Thus, the function of Dnmt1 of insects beyond its maintenance of cytosine methylation after DNA replication , and why its presence is variable, remains unclear.
One recently discovered role for the cytosine methylation machinery of insects is during gamete and embryo formation [3, 4, 8, 16, 22, 43, 48, 51]. In addition, for the large milkweed bug Oncopeltus fasciatus, Dmnt1 has been identified as critical during meiosis [4, 48]. This association with gamete formation may provide an insight into the variability of DNA methyltransferases of insects. Mechanisms of insect gamete formation is highly diverse [15, 20], and if Dnmt1 is related to meiosis during gametogenesis this may explain the diversity we see of this gene for insects. To test this idea, here we investigated if Dnmt1 is associated with meiotic pathways by investigating gene expression in a priori defined meiosis-related genes. We also use transcriptomics to compare Dnmt1 gene expression in knockdown vs control insects, and examine associations between Dnmt1 expression and other genes and pathways.
The large milkweed bug O. fasciatus has relatively high levels of cytosine methylation, as well as functional single copies of Dmnt1, Dmnt2, and Dmnt3 . Our previous studies investigated Dnmt1 during gamete formation of both ovaries and testes given these are rapidly dividing cells and Dnmt1 replicates cytosine methylation after a cellular division . When adult female O. fasciatus are injected with Dnmt1 RNA interference (RNAi), they fail to produce viable eggs and gamete production stops . When injected during a juvenile stage prior to formation of primary oocytes via meiosis in the developing ovary, Dnmt1 knockdown completely ablates gamete formation but the somatic ovaries develop normally . When Dnmt1 is knocked down during a pre-meiosis stage of juvenile testis development, resulting in a large reduction of cytosine methylation within the testes, males emerge as adults with significantly smaller testes that have fewer developing sperm . Similar to females, when injected as adults, spermatogenesis is blocked and no further sperm develop resulting in reduced fertility once sperm that had been in the process of development at the time of treatment have been used up . These results indicate that the reduction of Dnmt1 specifically affects production of gametes rather than gamete maturation. This led us to hypothesize that Dnmt1 influences gametogenesis at meiosis.
With this study, we used RNA-seq on RNAi and control injected O. fasciatus males to test the hypothesis that Dnmt1 influences spermatogenesis by influencing the expression of meiosis genes. In females, this reduction of fecundity is associated with the differential gene expression of several hundred genes . Like other insects, this differential gene expression is not associated with differential methylation [8, 13]. However, if Dnmt1 specifically affects the primary oocytes, which represent a small proportion of the cells within the adult ovary, it may be that any specific signatures of differential gene expression might be masked by the large number of unaffected somatic cells. The situation in the testis is very different. Each testis tubule is comprised of spermatocysts containing sperm at varying stages of development, and most cells within the testis are spermatogonia, spermatocytes, or developing spermatids. Given that the underlying phenotype of males and females appears to be similar, an impediment to the transition from germ cell (the diploid oogonia or spermatogonia) to gametes, we compared the gene expression pattern of testes following Dnmt1 knockdown using RNAi. For this species, targeting Dnmt1 with RNAi is very robust producing similar phenotypes and reductions of cytosine methylation for both males and females using multiple different constructs and controls [8, 48], including the injection timing and tissue used here. To capture any early events as opposed to more global changes that might result from downstream impacts of impeded development not directly related to the knockdown of Dnmt1, we sampled testes at 7-day and 14-day post-injection, which are earlier than the 21-day post-injection samples taken in the Washington et al. , where testes have significantly reduced numbers of developing spermatocytes and a highly disrupted structure. With these earlier sampling points, we aimed to capture the gene expression changes that were leading to the developmental arrest of spermatogenesis, presumably at the transition between diploid spermatogonia and haploid spermatocytes. We also tested our broad hypothesis of post-mitotic-specific arrest using a priori candidate genes representing different reproductive, molecular pathways and used all expressed gene to test if meiosis GO terms were enriched.
Dnmt1 RNAi knockdown
Dnmt1 was effectively knocked down in experimental males treated with ds-Dnmt1 RNA. Using qRT-PCR, we confirmed Dnmt1 was downregulated in the ds-Dnmt1 samples (7 days: mean fold expression difference = -0.677, t26 = 7.250, P = 8.609e-8; 14 days: mean fold expression difference = -0.679, t26 = 4.288, P = 0.00011). For both data sets (7 days, 14 days) Dnmt1 was statistically significantly downregulated in Dnmt1 RNAi treatment (Table 1).
Microscopy—Effect of Dnmt1 RNAi knockdown
Knockdown of Dnmt1 had a measurable phenotype in males sampled at both 7- and 14-day post-injection. Knockdown of Dnmt1 reduced the number of actively dividing spermatocysts (Fig. 1A). The number of spermatocysts that were positively stained with anti-pHH3 antibody, which indicates cell undergoing mitotic or meiotic division, was reduced in males at both 7-day post-knockdown (mean difference = − 3.45, t38 = 4.098, P = 0.0001; Fig. 1B) and 14-day post-knockdown (mean difference = − 2.1, t38 = 2.978, P = 0.0025; Fig. 1B).
a priori candidate genes
There were no differentially expressed a priori candidate genes at 7 days (Table 1). Six of seventeen were at 14 days (Table 1). However, they did not concentrate around a particular process. Three of six cell cycle control genes were differentially expressed, while two of five meiotic genes were differentially expressed.
a priori GO term enrichment
Using every expressed gene for analysis, none of the GO terms were enriched at 7 or 14 days after Dnmt1 knockdown (Table 2). In fact, contrary to our expectation, all of the five GO terms that we made predictions for were absent from the top of the list of genes rank-ordered by difference of expression at 7 days (i.e., they were statistically preferentially found at the bottom of the list of genes rank-ordered by their expression differences), while none were overrepresented at either the top or bottom of the rank-ordered list at 14 days (Table 2). The necrotic and apoptotic GO terms were not enriched at either sampling point.
Differential gene expression
We investigated gene expression at two sampling points—7 and 14 days—post-injection to identify the genes that are differentially expressed after knockdown of Dnmt1, which leads to the systematic depletion of spermatocysts among adult O. fasciatus testes (number of biological replicates: 13 7-day eGFP, 10 7-day dnmt1, 15 14-day eGFP, and 14 14-day Dnmt1). The individual treatments clustered together well at both sampling points (Fig. 2). The treatments are not differentiated at 7 days (Fig. 2A; Additional file 3), are but highly differentiated at 14 days (Fig. 2C; Additional file 4). There were 74 genes that were differentially expressed between the control treatment and the ds-Dnmt1 treated samples at 7-day post-injection (23 up-regulated; 51 down-regulated; Fig. 2B). There were 6,746 genes that were differentially expressed between the control treatment and the ds-Dnmt1 treated samples at 14-day post-injection (2,761 up-regulated; 3,985 down-regulated; Fig. 2D).
The overlap between the differentially expressed genes found here and those of Bewick et al.  using the same treatment against ovaries, was four for the 7-day list and ninety-six for the 14-day list. Both overlaps are statistically significant. The 7-day comparison had more than expected (mean simulated overlap = 1.32 P = 0.041) and the 14-day comparison had fewer than expected (mean simulated overlap = 119.24, P = 0.0022).
Gene Ontology (GO) term overrepresentation
Among the differentially expressed genes between Dnmt1 knockdown and eGFP control at 7-day samples, 113 GO terms were statistically significantly overrepresented. The top three Biological Processes terms with more than one gene were purine nucleobase biosynthetic process—GO:0009113, skeletal muscle organ development—GO:0060538, and mitotic sister chromatid cohesion—GO:GO:0007064 (P = 0.0014, 0.0016, 0.006, respectively); Molecular Function terms were calcium-dependent phospholipid binding—GO:0005544, cation binding—GO:0043169, and ribose phosphate diphosphokinase activity (P = 0.0017, 0.008, 0.0083, respectively); and for Cellular Compartments were nuclear meiotic cohesin complex—GO:0034991, Parkin-FBXW7-Cul1 ubiquitin ligase complex—GO:1,990,452, and integral component of Golgi membrane—GO:0030173 (P = 0.0087, 0.0087, 0.0125, respectively). Very few GO terms collapsed under a parent GO term in our semantic similarity analysis trying to understand what broad process was most perturbed with Dnmt1 knockdown. A complete term list can be found in Additional file 5.
For the 14-day samples, 322 GO terms were statistically significantly overrepresented. The top three Biological Processes terms with more than one gene were response to unfolded protein—GO:006986, double-strand break repair via nonhomologous end joining—GO:006303, and mRNA splice site selection—GO:0006376 (P = 0.00011, 0.00016, 0.00019, respectively), Molecular Function terms were nucleoside-triphosphatase activity—GO:0017111, single-stranded DNA binding—GO:0003697, and electron transfer activity—GO:0009055 (P = 0.00014, 0.00038, 0.00062, respectively), and for Cellular Compartments were ribonucleoprotein granule—GO:0035770, CCR4-NOT core complex—GO:0030015, sperm flagellum—GO:0036126 (P = 0.00015, 0.00025, 0.00025, respectively). Very few GO terms collapsed under a parent GO term in our semantic similarity analysis trying to understand what broad process was most perturbed with Dnmt1 knockdown. A complete term list can be found in Additional file 6.
For the overlap between the differentially expressed genes found here and those of Bewick et al.  for the 14-day samples, we found 168 GO terms. The top three Biological Processes terms with more than one gene were protein K63-linked deubiquitination—GO:0070536, melanocyte differentiation—GO:0030318, and negative phototaxis (P = 0.0002, 0.00025, 0.00032, respectively), Molecular Function GO terms were proteasome binding—GO:0070628, alditol:NADP + 1-oxidoreductase activity—GO:0004032, and estradiol 17-beta-dehydrogenase activity—GO:0004303 (P = 0.00012, 0.00029, 0.000221, respectively), and for Cellular Compartments were molybdopterin synthase complex—GO:0042629, apical cortex—GO:0045179, and outer dynein arm—GO:0036157 (P = 0.00017, 0.00691, 0.00072, respectively). Very few GO terms collapsed under a parent GO term in our semantic similarity analysis. A complete term list can be found in Additional file 7.
Producing gametes is achieved by a surprisingly diverse set of species- and lineage-specific set of mechanisms and we suggested that Dnmt1 may be part of this diversity. Our previous work with O. fasciatus shows that the genes for cytosine methylation are needed for successful gamete production [4, 8, 48] but the mechanistic underpinnings of this effect are unclear. Here, we reduced the gene expression of Dnmt1 of sexually mature O. fasciatus males. Given the phenotype observed here and previous studies, we predicted meiosis, and its molecular pathways would be highly perturbed, giving rise to a loss of meiotic progression and a cessation of spermatocyte production observed. Yet our candidate gene screen looking for differential expression of a priori predicted genes did not reveal an obvious molecular pathway. Both the collection of cell cycle and meiotic genes contained differentially expressed genes, but it was not the majority of either group; 3/6 for cell cycle and 2/5 for meiotic genes. These results align with the observation that cells appear to arrest at the meiotic stage, do not progress through further cell division which would require completed cell cycles, and that Dnmt1 might have a cytosine methylation independent function during reproduction.
As with the candidate gene analysis, our hypothesis received little support from our GO term enrichment analysis. None of the seven GO terms that we tested were enriched at 7 days (i.e., the genes with those GO terms were not preferentially found at the beginning of the rank-ordered gene list). In fact, the GO term searched were preferentially found at the end of the gene list at the 7-day mark, except for necrotic and apoptotic processes which did not have statistically significantly pattern (i.e., were among genes that showed the smallest differences of expression between the control and Dnmt1 knockdown samples). Like our candidate gene results, this suggests that few pathways were perturbed at this sampling point and there is no evidence that the cells are aborting themselves. At 14 days, none of the GO terms were enriched nor were they preferentially found at the end of the rank-ordered gene list. This means that these genes collectively increased their rank within the rank-ordered gene lists based on expression differences from 7 days to 14 days, but not at a rate that made them statistically enriched at the top of the 14-day gene list. The necrotic and apoptotic GO terms were not enriched at either timepoint suggesting that cells were not actively promoting either of these processes at either sampling point. Taken together, we suggest this supports the hypothesis that the gene expression difference at 14 days is a global response to Dnmt1 knockdown and is not highly targeted to any specific set of molecular pathways. In general, the cells do not have a detectable signal that they are undergoing any process to actively remove themselves from the cell population.
Given the lack of support for our hypothesis of a role for Dmnt1 in a meiosis pathway, we used the transcriptome to generate post hoc hypotheses. We found few gene expression difference at 7 days, but a large amount at 14 days post-treatment (~ 45% of expressed genes are differentially expressed). There was no clear consensus among the overrepresented GO terms to suggest how this cellular arrest is being achieved. All of the top five GO terms produced by the semantic similarity reduction for any three of the categories for both days were single standalone terms. The overlap between Bewick et al.  and our samples also did not produce any strong candidate pathways for Dnmt1’s influence after reduction. These pieces of evidence suggest that gametes are proceeding along a developmental pathway until Dnmt1 is needed, and when Dnmt1 is not present due to gene knockdown, the cell has a massive, pervasive, and non-specific disruption of its transcription environment. It appears that the differential gene expression seen is being driven by a healthy vs arrested cell contrast and that the effect of a reduction of Dnmt1 expression is a general, global response tied to chromatin condensation, not a set of specific genomic loci that are being targeted. Our RNA-seq results suggest that testes with reduced Dnmt1 expression have an increasingly degradation of their transcriptional environment. The 74 differentially expressed genes of the 7-day samples did not contain many transcription factors or developmental pathway genes as we expected. At 14 days, nearly half of all the expressed genes are differentially expressed—6794 of 14,984 genes (45.3%). Within ovaries treated in the same way and sampled at 10-day post-treatment, there were an intermediate amount of 236 differentially expressed genes . There was little overlap from our 7 days differentially expressed gene list with the ovary samples and fewer than statistically expected with the 14-day list.
We propose that Dnmt1’s action during gametogenesis influences chromosomal dynamics and nuclei condensation, based on multiple and repeatedly seen lines of evidence in our experiments. Gene expression knockdown of Dnmt1 does lead to reduced cytosine methylation [4, 8, 48]. However, cytosine methylation differences are not directly causal for differential gene expression for O. fasciatus , or for insects generally . This suggests that Dnmt1 actions to produce this phenotype are not driven by differences of cytosine methylation. When gene expression of Dnmt1 is reduced hundreds to thousands of gene are differentially expressed for both male and female reproductive tissues , this study). We posit that the differential gene expression seen between Dnmt1 gene expression knockdown and control samples is attributable to comparing samples, where one set is naturally progressing through gamete formation, while the other is in a highly arrested state. It does not appear that this arrest is an active process that the cells are undergoing, but rather is a function of the cell not being able to progress. This aligns with there being no signal from the GO term enrichment the cells are undergoing any cellular death process; necrotic or apoptotic; and that the cell looks healthy at a gross anatomical level. Even though differences of methylation of individual genomic loci do not lead to differences of gene expression between treatments, high cytosine methylation is associated with high gene expression between different loci [17, 18, 28, 49]. This high methylation–high expression association is particular to insects in contrast to plants, fungi, and vertebrates for which high cytosine methylation is often associated with heterochromatin . Highly condensed nuclei are seen for testis tubules with gene expression knockdown of Dnmt1 . This effect is most pronounced at the region, where spermatocysts transition from spermatogonia to spermatocytes through meiosis . When Dnmt1 gene expression is reduced by RNAi, it also leads to a reduction of cytosine methylation of somatic tissues . However, this reduction does not impact life span, suggesting that the lack of cytosine methylation within somatic tissues does not have a strong impact on the somatic cells of the organism . This also points to a difference between somatic and reproductive cells which aligns with the presence and absence of meiosis. What remains to be directly explained by this model is why reduced Dnmt1 gene expression knockdown or a lack of cytosine methylation is associated with highly condensed chromatin (or conversely an inability to de-condense chromatin) as gametogenesis progresses from earlier to later stages.
Careful molecular natural history has exposed much variation across insects for the level of gene body methylation [7, 39] and that differences of methylation are not casually associated with differences of gene expression . A role for Dmnt1 in chromosome dynamics addresses why its actions appear crucial yet independent function during reproduction, why cytosine methylation has persisted for insects, and why it is rarely associated with gene expression across the taxon. Even with this suggestion, it remains possible that cytosine methylation could directly influence gene expression in a subset of cells even if it does not systemically. This mechanism of action would explain why a reduction of Dnmt1 expression (1) leads to a highly perturbed transcriptional environment but why cytosine methylation is not associated with gene expression directly, (2) why Dnmt1 RNAi tissues have high condensed nuclei at the boundary, where spermatocysts usually enter the first meiotic division, and (3) why Dnmt1 has a specific and critical role during meiosis that likely is not driven by difference of DNA methylation between treatments.
To understand the functional role that Dnmt1 plays during gamete formation of male O. faciatus, we used RNAi to knockdown Dnmt1 gene expression at two sampling points post-treatment: 7 days and 14 days. After the treatment at two sampling points, we used confocal microscopy to phenotype the testes and determine the state of spermatocytes. We then profiled gene expression using RNA-seq of testes at the same two sampling points to understand how Dnmt1 knockdown perturbed the testes transcriptome.
Animal colony and husbandry
Oncopeltus fasciatus colony stocks were purchased from Carolina Biologicals (Burlington, NC). The colonies were maintained in Percival incubators under a 16:8 h light:dark cycle at 26 °C. The animals were fed organic raw sunflower seeds and had ad libitum deionized water. We needed to collect animals of known age so we removed eggs from the colonies and housed them in plastic storage containers with food and water. Nymphs were sexed at the fourth instar and separated into single sex colonies. Containers were checked daily for adult eclosions. New adults were placed into individual petri dishes with food and water.
RNA interference (RNAi) synthesis, administration, and quality control
We created RNAi constructs of Dnmt1 following Bewick et al. . Briefly, DNA templates of Dnmt1 were prepared via PCR. Following that, double stranded RNA was synthesized and then digested with an Ambion MEGAscript kit (ThermoFisher Sci, Waltham, MA) following the manufacturer’s protocol. This reaction was purified with a phenol:chloroform:IAA extraction followed by a sodium acetate precipitation. Concentration was measure by a Qubit using the ssRNA kit following the manufacturer’s protocol. Sense and anti-sense strands were then allowed to anneal. We used eGFP as an exogenous control construct. Before administration, the RNAi construct concentration was adjusted to 4 μg/μL in injection buffer (5 mM KCl, 0.1 mM NaH2PO4) for both Dnmt1 and eGFP. Sexually mature, virgin males (~ 7-day post-adult ecolsion) were injected with 3 µL ds-Dnmt1 RNA or eGFP using a pulled capillary needle between the third and fourth abdominal segments . Males were paired with a virgin female to allow for mating and stimulate sperm production. Pairs were kept in petri dishes under standard rearing conditions until they sampled. Males were haphazardly allocated to treatment group.
We have established previously there is no difference between controls using RNA constructs with no specificity to O. fasciatus genome sequence (e.g., eGFP) and buffer alone injections . Therefore, we did not include a buffer alone control. Previously, we have targeted RNAi constructs to two different regions of Dnmt1; the cytosine-specific DNA methyltransferase replication foci domain (RFD) and the DNA methylase domain (AdoMet). These give rise to identical phenotypes, which suggests minimal off-target effects . Thus, here, we only targeted against the RFD consensus domain.
RNAi validation by quantitative RT-PCR (qRT-PCR)
We assessed the effectiveness of Dnmt1 knockdown, using qRT-PCR. We synthesized cDNA using 500 ng RNA with qScript cDNA SuperMix (Quanta Biosciences, Gaithersburg, MD). Validated primers were used from Bewick et al. . GAPDH was used as the endogenous reference gene. We used Roche LightCycler 480 SYBR Green Master Mix with a Roche LightCycler 480 (Roche Applied Science, Indianapolis, IN) with 3 technical replicates of 10 μL reactions as previously reported . We used the ΔΔCT method to estimate differences of expression using eGFP samples as our comparison group . We had 14 7-day eGFP, 14 7-day Dnmt1, 14 14-day eGFP, and 13 14-day Dnmt1 biological replicates.
Phenotypic analysis of Dnmt1 knockdown
Virgin, adult males were collected on the day of emergence. Males were injected with ds-eGFP or ds-Dnmt1 as described above at 7–10-day post-adult emergence. Males were haphazardly assigned to two dissection days. One group of males were dissected 7-day post-RNAi treatment and a second group was dissected 14-day post-RNAi treatment.
Testes were dissected from the males in PBS and processed for microscopy as described below. To assess the activity of early stage spermatocytes, the testis tubules were stained with anti-phosphohistone H3 Ser 10 antibody (pHH3; Millipore antibody 06-570, Sigma Aldrich, St. Louis, MO). Anti-pHH3 stains for chromosome condensation in preparation for mitosis and meiosis [19, 38]. Primary antibody staining was visualized with a secondary antibody, Alexa Fluor goat–anti-rabbit 647 (ThermoFisher Scientific, Waltham, MA). Testis tubules were counterstained with DAPI (0.5 μg/mL PBT) to visualize nuclei. Negative controls in which primary antibody was absent showed no non-specific binding of the secondary antibody (data not shown.) Testis tubules were visualized on a Zeiss 710 confocal microscope or an EVOS (ThermoFisher Scientific, Waltham, MA).
To quantify the number of actively dividing spermatocysts, testis tubules from males were examined. Anti-pHH3 antibody stains both spermatogonia undergoing mitotic divisions and spermatocysts undergoing meiotic divisions (Fig. 1A). The numbers of positively stained spermatocysts in either control (ds-eGFP) or knockdown (ds-Dnmt1) males were compared at within each sample using a t test with unequal variance and the prediction that ds-Dnmt1 samples would have fewer dividing cells (n = 20 per treatment, except 7-day Dnmt1 had 18).
The testes of males of each treatment were dissected out in ice-cold 1X PBS at the appropriate sampling point, flash frozen with liquid nitrogen, and stored at − 80 °C until nucleotide extraction. Total RNA was extracted using a Qiagen Allprep DNA/RNA Mini Kit (Qiagen, Venlo, The Netherlands) following the manufacturer’s protocol. Homogenization of the testes were performed with a handheld Kimble pellet pestle in RLT buffer. Quantification was done with a Qubit fluorometer using the RNA BR kit.
RNA-seq high-throughput library preparation and sequencing
The extracted RNA of O. fasciatus ds-Dnmt1 and control (ds-eGFP) biological replicates at 7- and 14-day post-injection was used to construct poly-A selected Illumina TruSeq Stranded RNA LT Kit (Illumina, San Diego, CA) following the manufacturer’s instructions with limited modifications. The starting quantity of total RNA was adjusted to 1.3 µg, and all reagent volumes were reduced to a third of the described quantity. We targeted 10 M 2 × 150 bp read pairs per biological replicate using an Illumina NextSeq 500 with v3.1 chemistry. There were 13 7-day eGFP, 11 7-day Dnmt1, 16 14-day eGFP, and 15 14-day Dnmt1 libraries. Libraries were sequenced at the Georgia Genomics and Bioinformatics Core (Athens, GA, USA).
Read quality control and mapping
Reads were initially assessed for quality with fastQC (v0.11.9; default settings;  to establish a baseline. Reads had adapters trimmed with cutadapt (v2.8,–trim-n -O 3 -u -2 -U -2 -q 10,10 -m 30;  using the TruSeq adapter sequences. Reads were again assessed with fastQC with default settings. Overlapping reads were combined with FLASh (v 1.2.11, default settings; . As a final QC step, reads that mapped to rRNA genes were removed with SortMeRNA (v4.3.3, all gff entries annotated as rRNA) .
We used the NAL i5k O. faciatus genome (the “BCM-After-Atlas” version) and annotation (Official Gene Set v1.2; . These were downloaded from the NAL i5k site: https://i5k.nal.usda.gov/content/data-downloads. This was the most current version of the genome and gene annotation at the time of analysis. HISAT2 (v2.2.1; no soft-clipping;  was used to map reads to the genome (Additional file 1. Extended reads (i.e., reads that were combined by FLASh had ~ 30% higher mapping rate. Mappings were converted to read counts by StringTie (v2.1.7;  following the manual instructions for export to DESeq2.
We updated the functional annotation of the O. faciatus proteome/transcriptome using eggNOG-mapper webserver at the level of Insecta (http://eggnog-mapper.embl.de/; v 2.19; . This annotated 12,526 of 19,615 gene models (63.8%) with Gene Ontology terms, which allowed us to perform a GO term enrichment analysis and a GO term overrepresentation analysis.
Differential gene expression
Read counts were imported into R using tximport (Bioconductor v1.20.0;  following the manual’s instructions. We used R (v4.1.0 , within an RStudio IDE (build 492 , for the differential gene expression analysis.
DESeq2 (Bioconductor v1.32.0; default settings;  was used for all DGE analyses following the programmers’ suggestions for exploratory data analysis and sample/analysis quality control. eGFP was set and used as the comparison group for all analyses. 7-day and 14-day samples were analyzed separately as previously discussed to generate a list of differentially expressed genes between control and experimental treatments. The model matrix specification was also checked to ensure correct specification (i.e., that program was contrasting the samples in the correct way). After importing and initial analysis, samples were plotted with a PCA to visually check for outliers according to the manual’s recommendation. One 7-day Dnmt1 sample was removed. Two 14-day samples were also removed—one Dnmt1 and one eGFP. All were removed due to poor library quality.
After removal of the outliers each analysis was repeated using the same settings of DESeq2 and the results were again quality controlled. We used the default dispersion estimator and shrinkage method, apeglm . We used s values to estimate statistical significance after false discovery rate correction at the level of 0.05 . Results were visualized with the fviz_pca_ind function of the factoextra R package .
We also compared the overlap of the differentially expressed genes here to the differentially expressed genes of O. fasciatus ovaries under the same treatment scheme using the intersect function of R. A simulation analysis was conducted to see if any overlap observed was statistically significantly enriched or depleted .
a priori Candidate gene screen
The first test of our hypothesis was to assess the influence Dnmt1 knockdown had on a series of candidate genes. These genes were chosen based on literature searches for genes involved in meiosis and spermatogenesis of insects. These include cell fate and proliferation genes that we did not expect to be differentially expressed (members of the Wnt and Frizzled families), cell cycle control genes some of which we did expect to be differentially expressed (Vasa) and some of which did not expect to be differentially expressed (Cyclin B3, Cyclin D2, Cyclin D3, Cell Division Cycle 20, Cell Division Cycle 25 family members), maintenance of chromosome genes that we did expect to be differentially expressed (Structural Maintenance of Chromosome 3 family members), meiotic transition gene that we did expect to be differentially expressed (Boule), and meiotic recombination genes that we did expect to be differentially expressed (SPO11 Initiator of Meiotic Double Stranded Breaks, MutS Homolog 5, Meiotic Nuclear Divisions 1, Homologous-Pairing Protein 2). None of the candidates had specific directional changes, except decreased expression of Vasa and Boule within the Dnmt1 knockdowns. We extracted the raw P values for expression differences from the DESeq2 results and compared them with FDR corrected P values we generated using the Benjamini–Hockberg procedure .
a priori GO term enrichment analysis
As secondary test of our hypothesis, we selected five high-level GO terms that represented the cellular processes that we thought were being perturbed or not after Dnmt1 knockdown. These were GO:0007049 Cell Cycle, GO:0000278 Mitotic Cell Cycle, GO:0051321 Meiotic Cell Cycle, GO:0007276 Gametogenesis, and GO:0007283 Spermatogenesis. With our previous observed phenotypes and the one observed with the microscopy here, we expect Meiotic Cell Cycle, Gametogenesis, and Spermatogenesis to be enriched. We expected Cell Cycle and Mitotic Cell Cycle not to show enrichment. We included GO:0070265 Necrotic Cell Death and GO:0006915 Apoptotic Process to assess if cells were undergoing these processes or if they had simply arrested. We used S values from DESeq2 to rank order all expressed genes at 7 and 14 days separately. We used fgsea R package  with default parameters, except allowing 1,500 genes as a maximum, to test if genes annotated with these GO terms were enriched at the top of these lists. Results were visualized with the same R package.
Gene ontology (GO) term overrepresentation analysis
We found overrepresented terms to give us biological processes that might be perturbed using the topGO package of R (v2.48.0; . This analysis is often termed an enrichment analysis , but statistically it tests if there is an overrepresentation of a GO term from a list of interesting genes (usually those that are differentially expressed compared with the expectation of the same number of random picks. It does not calculate if a specific GO term (or another type of annotation, e.g., KEGG is enriched at the beginning of an ordered gene set. We performed GO term tests using Fisher’s exact test with the weighted algorithm. GO terms of genes of interest were compared to all expressed genes within our samples. Updated GO terms can be found in Additional file 2.
We performed a semantic similarity analysis using REVIGO using default settings  to refine our overrepresented GO terms to even higher level summaries of the processes involved and identify central processes enriched from our treatments.
Availability of data and materials
All high-throughput data are available under NCBI BioProject # PRJNA957711. Data for testes tube nuclei are provided in Additional file 8.
Agrelius TC, Altman J, Dudycha JL. The maternal effects of dietary restriction on Dnmt expression and reproduction in two clones of Daphnia pulex. Heredity. 2023;130:73–81.
Alexa A, Rahnenfuhrer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006;22:1600–7.
Arsala D, Wu X, Yi SV, Lynch JA. Dnmt1a is essential for gene body methylation and the regulation of the zygotic genome in a wasp. PLoS Genetics. 2022;18:e1010181.
Amukamara AU, Washington JT, Sanchez Z, McKinney EC, Moore AJ, Schmitz RJ, et al. More than DNA methylation: does pleiotropy drive the complex pattern of evolution of Dnmt1? Frontiers Ecol Evol. 2020;8(10):3389.
Andrews S (2010) FastQC: A quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Benjamini Y, Hockberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995;57:289–300.
Bewick AJ, Vogel KJ, Moore AJ, Schimitz RJ. Evolutin of DNA methylation across insects. Molec Biol Evol. 2017;34:654–65.
Bewick AJ, Sanchez Z, McKinney EC, Moore AJ, Moore PJ, Schmitz RJ. Dnmt1 is essential for egg production and embryo viability in the large milkweed bug Oncopeltus fasciatus. Epigenetics Chromatin. 2019;12:1–14.
Cantalapiedra CP, Hernández-Plaza A, Letunic I, Bork P, Huerta-Cepas J. eggNOG-mapper v2: functional annotation, orthology assignments, and domain prediction at the metagenomic scale. Molec Biol Evol. 2021;38:5825–9.
Chesebro J, Hrycaj S, Mahfooz N, Popadic A. Diverging functions of scr between embryonic and postembryonic development in a hemimetabolous insect Oncopeltus fasciatus. Develop Biol. 2009;329:142–51.
Cunningham CB, Douthit MK, Moore AJ. Octopaminergic gene expression and flexible social behaviour in the subsocial burying beetle Nicrophorus vespilloides. Insect Mol Biol. 2014;23:391–404.
Cunningham CB, Ji L, McKinney EC, Benowitz KM, Schmitz RJ, Moore AJ. Changes of gene expression but not cytosine methylation are associated with male parental care reflecting behavioural state, social context and individual flexibility. J Exp Biol. 2019;222:jeb188649.
Duncan EJ, Cunningham CB, Dearden PK. Phenotypic plasticity: what has DNA methylation got to do with it? Insects. 2022;13:110.
Engelhardt J, Scheer O, Stadler PF, Prohaska SJ. Evolution of DNA methylation across Ecdysozoa. J Mol Evol. 2022;90:56–72.
Ewen-Campen B, Schwager EE, Extavour CGM. The molecular machinery of germ line specification. Mol Reprod Develop. 2010;77:3–18.
Gegner J, Gegner T, Vogel H, Vilcinskas A. Silencing of the DNA methyltransferase 1 associated protein 1 (DMAP1) gene in the invasive ladybird Harmonia axyridis implies a role of the DNA methyltransferase 1-DMAP1 complex in female fecundity. Insect Molec Biol. 2020;29:148–59.
Gladstad KM, Hunt BG, Goodisman MAD. Evolutionary insights into DNA methylation in insects. Cur Opin Insect Sci. 2014;1:25–30.
Glastad KM, Goodisman MA, Yi SV, Hunt BG. Effects of DNA methylation and chromatin state on rates of molecular evolution in insects. G3: Genes Genomes Genetics. 2016;6:357–63.
Hans F, Dimitrov S. Histone H3 phosphorylation and cell division. Oncogene. 2001;20:3021–7.
Hemming BS. Insect development and evolution. New York: Cornell University Press; 2018.
Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13.
Ivasyk I, Olivos-Cisneros L, Valdés-Rodríguez S, Doual M, Jang H, Schmitz RJ, Kronauer DJC. DNMT1 mutant ants develop normally but have disrupted oogenesis. Nat Commun. 2023;14:2201.
Kassambara A, Mundt F (2020) Factoextra: Extract and visualize the results of multivariate data analyses. v1.0.7. CRAN.R-project.org/package=factoextra
Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotech. 2019;37:907–15.
Kopylova E, Noé TH. SortMeRNA: Fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211–7.
Korotkevich F, Sukhov V, Budin N, Shpak B, Artyomov MN, Sergushichev A. Fast gene set enrichment analysis. BioRxiv. 2021. https://doi.org/10.1101/060012.
Law JA, Jacobsen SE. Establishing, maintaining, and modifying DNA methylation patterns in plants and animals. Nature Rev Genetics. 2010;11:204–20.
Libbrecht R, Oxley PR, Keller L, Kronauer DJC. Robust DNA methylation in the clonal raider ant brain. Curr Biol. 2016;26:391–5.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25:402–8.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biolo. 2014;15:550.
Magoc T, Salzberg S. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics. 2011;27:2957–63.
Maleszka R, Kucharski R. Without mechanism, theories and models in insect epigenetics remain a black box. Trends Genet. 2022;38:1108–11.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet Journal. 2011;17:10–2.
Morandin C, Brendel VP. Tools and applications for integrative analysis of DNA methylation in social insects. Mole Ecol Res. 2021;22:1656–74.
Panfilio KA, Jentzsch IMV, Benoit JB, Erezyilmaz D, Suzuki Y, Colella S, et al. Molecular evolutionary trends and feeding ecology diversification in the Hemiptera, anchored by the milkweed bug genome. Genome Biol. 2019;20:64.
Patalano S, Alsina A, Gregorio-Rodriguez C, Bachman M, Dreier S, Hernado-Herraez I, Nana P, Balasubramanian S, Sumner S, Reik RS. Self-organization of plasticity and specialization in a primitively social insect. Cell Syst. 2022;13:768–79.
Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg S. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature Biotechnol. 2015;33:290–5.
Prigent C, Dimitrov S. Phosphorylation of serine 10 in histone H3, what for? J Cell Sci. 2003;116:3677–85.
Provataris P, Meusemann K, Neihuis O, Garth S, Misof B. Signatures of DNA Methylation across insects suggest reduced DNA methylation levels in holometabola. Genome Biol Evol. 2018;10:1185–97.
R Core Team (2021) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
RStudio Team (2022) RStudio: Integrated Development for R. RStudio, PBC, Boston, MA URL http://www.rstudio.com/.
Schmitz RJ, Lewis ZA, Goll MG. DNA methylation: shared divergent features across eukaryotes. Trends Genet. 2019;35:818–27.
Schulz NKE, Wagner CI, Ebeling J, Raddatz G, Diddens-de Buhr M, Lyko F, et al. Dnmt1 has an essential function despite the absence of CpG DNA methylation in the red flour beetle Tribolium castaneum. Sci Reports. 2018;8:16462.
Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2015;4:1521.
Stephens M. False discovery rates: a new deal. Biostatistics. 2016;18:2.
Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE. 2011;6: e21800.
Vogt G. Evolution, functions, and dynamics of epigenetic mechanisms of annimals. In handbook of epigenetics. Cambridge: Academic Press; 2022.
Washington JT, Cavender KR, Amukamara AU, McKinney EC, Schmitz RJ, Moore PJ. The essential role of Dnmt1 in gametogenesis in the large milkweed bug Oncopeltus fasciatus. Elife. 2021;10:62202.
Yan H, Bonasio R, Simola DF, Liebig J, Berger SL, Reinberg D. DNA methylation in social insects: how epigenetics can control behavior and longevity. Annu Rev Entomol. 2014;60:435–52.
Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics. 2018;35:2084–92.
Zwier MV, Verhulst EC, Zwahlen RD, Beukeboom LW, van de Zande L. DNA methylation plays a crucial role during early Nasonia development. Insect Molec Biol. 2012;21:129–38.
We appreciate the discussions and contributions from other members of our laboratory, especially Ahva Potticary.
This work was funded by the USDA–ARS Non-Assistance Cooperative Agreement “Managing whiteflies and whitefly transmitted viruses in vegetable crops in the southeastern U.S.” (#58-6080-9-006) to AJM and support from the National Science Foundation (MCB-1856143) to RJS.
Ethics approval and consent to participate
The authors have no conflicts of interest to declare.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Contains the metadata for the RNA-seq libraries used here.
Contains the updated GO term annotations for Oncopeltus fasciatus.
Contains the output from DESeq2 for the differential gene expression of the 7-day samples.
Contains the output from DESeq2 for the differential gene expression of the 14-day samples.
Contains the output from topGO for the GO term overrepresentation analysis of the 7-day samples.
Contains the output from topGO for the GO term overrepresentation analysis of the 14-day samples.
Contains the output from topGO for the GO term overrepresentation analysis of the Bewick et al.  and 14-day samples overlap
Contains the raw data for figure one, nuclei count of dividing gametes.
About this article
Cite this article
Cunningham, C.B., Shelby, E.A., McKinney, E.C. et al. The role of Dmnt1 during spermatogenesis of the insect Oncopeltus fasciatus. Epigenetics & Chromatin 16, 28 (2023). https://doi.org/10.1186/s13072-023-00496-5