Overview

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. [8]. 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 [10]. 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 [8]. 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 [8]. 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. [8]. 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 [11]. We used the ΔΔCT method to estimate differences of expression using eGFP samples as our comparison group [29]. 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).

RNA extraction

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; [5] 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; [33] 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; [31]. As a final QC step, reads that mapped to rRNA genes were removed with SortMeRNA (v4.3.3, all gff entries annotated as rRNA) [25].

We used the NAL i5k O. faciatus genome (the “BCM-After-Atlas” version) and annotation (Official Gene Set v1.2; [35]. 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; [24] 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; [37] following the manual instructions for export to DESeq2.

Functional annotation

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; [9]. 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; [44] following the manual’s instructions. We used R (v4.1.0 [40], within an RStudio IDE (build 492 [41], for the differential gene expression analysis.

DESeq2 (Bioconductor v1.32.0; default settings; [30] 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 [50]. We used s values to estimate statistical significance after false discovery rate correction at the level of 0.05 [45]. Results were visualized with the fviz_pca_ind function of the factoextra R package [23].

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 [12].

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 [6].

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 [26] 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; [2]. This analysis is often termed an enrichment analysis [21], 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 [46] to refine our overrepresented GO terms to even higher level summaries of the processes involved and identify central processes enriched from our treatments.