Effect of heritable symbionts on maternally-derived embryo transcripts

Maternally-transmitted endosymbiotic bacteria are ubiquitous in insects. Among other influential phenotypes, many heritable symbionts of arthropods are notorious for manipulating host reproduction through one of four reproductive syndromes, which are generally exerted during early developmental stages of the host: male feminization; parthenogenesis induction; male killing; and cytoplasmic incompatibility (CI). Major advances have been achieved in understanding mechanisms and identifying symbiont factors involved in reproductive manipulation, particularly male killing and cytoplasmic incompatibility. Nonetheless, whether cytoplasmically-transmitted bacteria influence the maternally-loaded components of the egg or early embryo has not been examined. In the present study, we investigated whether heritable endosymbionts that cause different reproductive phenotypes in Drosophila melanogaster influence the mRNA transcriptome of early embryos. We used mRNA-seq to evaluate differential expression in Drosophila embryos lacking endosymbionts (control) to those harbouring the male-killing Spiroplasma poulsonii strain MSRO-Br, the CI-inducing Wolbachia strain wMel, or Spiroplasma poulsonii strain Hyd1; a strain that lacks a reproductive phenotype and is naturally associated with Drosophila hydei. We found no consistent evidence of influence of symbiont on mRNA composition of early embryos, suggesting that the reproductive manipulation mechanism does not involve alteration of maternally-loaded transcripts. In addition, we capitalized on several available mRNA-seq datasets derived from Spiroplasma-infected Drosophila melanogaster embryos, to search for signals of depurination of rRNA, consistent with the activity of Ribosome Inactivating Proteins (RIPs) encoded by Spiroplasma poulsonii. We found small but statistically significant signals of depurination of Drosophila rRNA in the Spiroplasma treatments (both strains), but not in the symbiont-free control or Wolbachia treatment, consistent with the action of RIPs. The depurination signal was slightly stronger in the treatment with the male-killing strain. This result supports a recent report that RIP-induced damage contributes to male embryo death.

Embryo collection. Approximately 40-50 three-day-old virgin females, from each replicate of each treatment, were allowed to mate in cages with Wolbachia-free (W−) CS males during the collection period, and allowed to lay eggs on cornmeal food plates. The initial batch of eggs was discarded to improve the chances of collecting fertilized eggs from the same stage 44 . Thereafter, egg laying was monitored and cornmeal plates were changed approximately every 45 min, so as to collect embryos that were on average ~60-75 min old. Eggs were collected from each replicate with a small brush 6-8 times over a 2-day period for each treatment and the control. The eggs were placed in sterile RNase-free 1.7 ml microtubes, and immediately put on dry ice during the collection period, after which they were transferred to −80 °C for storage. RNA extraction, library preparation, and sequencing. Three to four biological replicates per condition (MSRO, Hyd1, wMel, and control) were used for the extractions (see Fig. 1). RNA was extracted per collection tube of eggs (mentioned above) with Trizol ® Plus RNA Purification System (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol. RNA from tubes that belonged to the same biological replicate within a treatment was pooled. All RNA samples were DNase-treated with Ambion ® DNA-free (Invitrogen) to remove any DNA contamination. Total RNA was quantified with a NanoDrop ® ND-1000 UV spectrophotometer (NanoDrop Technologies, Wilmington, DE), and sample quality and integrity were further tested with the Agilent 2100 BioAnalyzer (Agilent Inc., Santa Clara, CA).
Total RNA was submitted to the Texas AgriLife Genomics and Bioinfomatics Services facility for library preparation, sequencing, and demultiplexing. Three biological replicates were subjected to multiplex NuGen  (4), whereas underlined samples were pooled into another lane (5). Library labelled "Sample_12_TruSeq_12" is the result of combining total RNA from Hyd1 biological replicates 2 and 3, and thus considered a technical, rather than biological, replicate. This library was excluded from the differential expression analyses.
www.nature.com/scientificreports www.nature.com/scientificreports/ library preparation and pooled into a single sequencing lane of the Illumina (San Diego, CA) Genome Analyzer II platform (single-end; 76 bp read length; italicized samples in Fig. 1). The remaining biological replicates were subjected to the Illumina mRNA TruSeq kit library preparation protocol (12 libraries

Analyses of differential expression (DE). Raw reads have been deposited in the NCBI SRA Database
under Accession Numbers SRR7279355-SRR7279369 (BioProject PRJNA474708; BioSample SAMN09370647). Command lines, as well as input and output files needed to replicate our analyses, from the different pipelines used are available in Supporting Data Files S1-S4. Raw RNA sequence files were first processed with Trimmomatic (v 0.35) 45 to remove adapters and low quality reads (see Fig. 1). The processed runs were then mapped to the Flybase Drosophila genome (version dmel-all-r6.18) with HiSat2 (v.2.0.5) 46 to obtain bam files. The files containing the mapped read information were then analysed with HTseq (0.6.1) 47 to obtain read counts for both genes and exon regions specified in a corresponding gff file for the dmel-all-r6.18 genome sequence. The resulting count files were then manually formatted into a count matrix suitable for differential expression analysis (Supporting Data File S2).
Differential expression analyses were performed with three different pipelines edgeR v.3.20.5 48 and Limma + Voom v. limma_3.34.5 49 and DeSeq2 v.1.18.1 50 . Command line and relevant input files for these analyses are provided in Supporting Data Files S1-S3. All programs conducted their respective tests on genes/exons with 10 or more mapped reads per replicate within at least one of two the treatments being compared. For the EdgeR pipeline both the glm and QFglm models were used for DE. For Deseq2, the wald statistical test was used for DE analysis. Before analysis in both Limma-Voom and EdgeR pipeline, two genes with extremely high counts, 16S rRNA (FBgn0013686) and ef-1α (FBgn0284245), were removed. Genes with an adjusted p-value (or false discovery rate; FDR) of 0.05 or less were considered to be potentially differentially expressed (DE). Power analyses for differential expression. To identify possible limitations of our experimental design in the detection of DE genes, we performed a statistical power analysis with the method of Hart et al. 51 . The parameters (and rationale) used to run the simulations are provided in Table 1 Table 1. Parameters used to run the simulations of the power analyses. All possible combinations of these parameters (60 total) were run. a based on the Lander/Waterman estimate of coverage for our data; i.e., mean coverage of 23X per replicate (range 19-28 among our replicates). b based on our minimum number of reads per replicate (in one treatment) required for a DE test (see text) c These values encompass the range of most of the values of Biological Coefficient of Variation (BCV; estimated by edgeR) obtained (i.e., a "Common" value of ~0.31 and a stabilized "Trend" value of ~0.2; see Fig. 2). d reflects the range of replicates per treatment used in this study e power was high and stable beyond an effect size of 4 (see Fig. 5). Thus, higher effect sizes are not reported. www.nature.com/scientificreports www.nature.com/scientificreports/

Analyses of depurination signal in the sarcin-ricin loop (SRL) of the 28 S rRNA of Drosophila. RIP
toxins remove a specific adenine present in the sarcin-ricin loop (SRL) of the 28 S rRNA leaving an abasic site (i.e., the backbone remains intact) 52 . When a reverse transcriptase encounters an abasic site, it preferentially adds an adenine in the nascent complementary (c)DNA strand 53 . This property, which results in an incorrect base at the RIP-depurinated site in the cDNA and all subsequent PCR amplification steps, can be used to detect evidence of RIP activity in any procedure that relies on reverse transcription (e.g. RNA-seq or reverse-transcription qPCR). To examine whether a signal of depurination consistent with RIP activity was detectable in Spiroplasma-infected flies, we used the RNA-seq data generated in the present study, as well as RNA-seq data (also derived from poly-A-tail-enriched RNA libraries) from two published studies that also compared patterns of gene expression from Spiroplasma-infected and Spiroplasma-free D. melanogaster embryos: NCBI Acc. No. PRJDB4469 30 ; and PRJNA318373 35 . The Spiroplasma strain examined by 35 54 (for the prokaryotic genomes) to attempt mapping the "28S reads" to the following reference genomes: (1) the D. melanogaster genome (BDGP6); (2) MSRO-UG (NZ_PETG00000000.1); (3) a draft genome of Spiroplasma Hyd1 (unpublished data); and (4) Wolbachia wMel (NC_002978). As expected, the "28S reads" only mapped to (1), and specifically to regions whose annotation terms included "28S rRNA". To visualize and count the shift from A to T (or other bases) in the source RNA pool, the "28S reads" were mapped again to the NR_133562.1:2591-3970 in Geneious v.11.1.2 (Biomatters Inc., Newark, NJ; "low sensitivity mode"; maximum gap size = 3; iterate up to 25 times). The number of reads containing each of the four bases or a gap at the target site was counted by selecting the position at all the reads to be counted, and recording the counts reported by Geneious under the "Nucleotide Statistics" option (gapped reads were excluded from subsequent analyses). The proportion of reads with an A at the target site (i.e., putatively intact rRNA) was calculated and compared among treatments and replicates. We used a mixed model in JMP Pro 13 (SAS Institute Inc., Cary, NC) to examine the effect of symbiont treatment (fixed) on the proportion of adenines (arcsine square root transformed) at the target position (excluding gaps from total number of reads). Source Study was treated as a random effect.

Results
The present study aimed at determining whether heritable symbionts alter the composition of maternally-loaded mRNAs in D. melanogaster. To achieve this, we used mRNA-seq to evaluate differential expression (DE) in Drosophila embryos lacking endosymbionts (control) to those harbouring one of the following heritable endosymbionts: the male-killing Spiroplasma poulsonii strain MSRO-Br; the CI-inducing Wolbachia strain wMel; and the non-male-killing Spiroplasma poulsonii strain Hyd1 (native to D. hydei). A power analysis was performed to determine the limitations of our experimental design. A secondary goal of this study capitalized on several available mRNA-seq datasets derived from Spiroplasma-infected Drosophila melanogaster embryos, to search for signals of damage (i.e., depurination) to Drosophila rRNA, consistent with the activity of Ribosome Inactivating Proteins (RIPs), which are encoded in the genomes of several Spiroplasma strains. This assay examined the proportion of adenines at the RIP target position (indicative of intact rRNA) in Spiroplasma-infected vs. Spiroplasma-free treatments. Differential expression. The number of reads that was obtained, passed QC, and mapped to the Drosophila genome is shown in Table 2. The PCA plot did not reveal any particular grouping of samples according to www.nature.com/scientificreports www.nature.com/scientificreports/ treatment (Fig. 3). PC1 explained 51% of the variance and revealed a slight "batch" effect, where the "GAII" replicates (i.e., 1_gall, 2_gall, and 3_gall) fell to the left of the other replicates. Nonetheless, removal of these replicates did not lead to a better grouping of replicates within each treatment (not shown). Thus, we considered that exclusion of the "GAII" replicates was not justified. Inclusion or exclusion of the Hyd1 technical replicate (e.g. n = 4 vs. 3 replicates, respectively) did not affect the differential expression analyses (only n = 3 is shown).
All differential expression (DE) tests were performed by comparing each symbiont treatment to the symbiont-free control. None of the eight pipelines used detected differentially expressed (DE) genes in the MSRO (male-killing Spiroplasma) and wMel (CI-inducing Wolbachia) treatments (Table 3). For the Hyd1 (non-male-killing Spiroplasma) treatment, half of the pipelines did not detect DE genes, whereas the other half detected a few DE genes ( Table 3). The Deseq2 pipeline and the edgeR pipeline (glm model only) detected 5-12 differentially expressed (DE) genes in the Hyd1 treatment based on 'by gene' and 'by exon' analyses (Table 3). Only five DE genes (ect, Osi6, Osi7, FBgn0038339, and FBgn0037099) were found in common among these pipelines. For all of these, the Hyd1 treatment exhibited lower expression, and lower variation among replicates, than the other treatments (Table 3 and Fig. 4). Fig. 5 and the rationale for the parameters assumed is provided in Table 1. Assuming a coverage of 20X and a coefficient of variation of 0.3, we should expect to detect: ~100% of genes differentially expressed by log-fold changes ≥3; at least 60% of genes with log-fold changes ≥2, and at least ~40% of genes with log-fold changes ≥1.75. Therefore, these results suggest that our design and data have sufficient power to detect log-fold changes ≥3, but more limited power below that log-fold change value.

Power analyses for differential expression. A plot of the power analyses is shown in
Signals of depurination. The number of reads mapped to the target site, as well as the proportion of adenines (representative of non-depurinated rRNA) by treatment replicate, are shown Fig. 6 and Table S1. Symbiont treatment had a significant effect on the proportion of adenines (F 4,26 = 29.8; p < 0.0001), whereas source study did not (Wald p = 0.4895). Both MSRO (UG and BR) treatments had a significantly (p ≤ 0.001) smaller proportion of As (mean = 96.78%) than the control (mean = 99.79%) and wMel (mean = 100%) treatments. Hyd1 exhibited intermediate proportions of As (mean = 98.53%; Fig. 6). In other words, we detected a mean depurination of 3.22% for MSRO and 1.47% for Hyd1, and effectively no depurination for the control or wMel treatment. Pooling of the two substrains of MSRO (i.e., MSRO-UG and MSRO-BR) into one treatment produced similar results (not shown). The Harumoto et al. 30 dataset, which included two embryo stages and differentiated between male and female embryos, did not exhibit substantially different depurination levels between stages or sex (e.g. mean %A was 98.27 for females and 98.49 for males; see Table S1).
Outcome of mapping mRNAseq reads to symbiont genomes. To determine the degree to which poly-A-tailed enriched RNAseq data yield reads assignable to the bacterial symbionts, we used Bowtie2 (v.2.3.4) 54 with default settings to map the reads from each symbiont-infected replicate to the closest symbiont reference genome available. The number of reads per replicate mapped to the symbiont genome ranged from 5 (for a Hyd1 replicate) to 8619 (for a wMel replicate; Table S2), representing a very small fraction of the total reads. The percentage of these putative symbiont-derived reads that mapped to the symbiont ribosomal RNA genes was quite variable (range 3.94-100%; mean 52%). The above numbers of symbiont-derived genes are inadequate for analyses of differential expression (DE) of the symbiont. Collectively these findings suggest that libraries generated via a poly-A-tail enrichment process tend to yield a very low proportion of reads of bacterial origin; a result that is expected given that polyadenylation is not part of the process of bacterial mRNA maturation.

Discussion
Effects of symbionts on early embryo mRNA composition. The present study aimed to test whether infection by a heritable bacterium strain that kills males during the embryonic stage influences the composition of mRNA transcripts present in the early embryo (ca. 60-75 min post-oviposition), and thus, before the onset of zygotic transcription. Based on several analysis pipelines and assumptions, we did not detect any differential expression in genes among embryos harbouring the male-killing Spiroplasma strain (MSRO-Br), the CI-inducing Wolbachia strain, and the symbiont-free control. These results suggest that these symbionts do not influence the composition of initial maternally-loaded transcripts or their degradation up to the examined stage, unless they do so to a degree below what was detectable by our experimental design (e.g. below a log-fold change of 3; see Power Analyses). Furthermore, it is possible, that such cytoplasmically transmitted symbionts could influence the composition of maternally-loaded proteins, or regulate the protein complement of the early embryo in a transcriptionally-independent manner (e.g. through post-transcriptional or post-translational controls). To our knowledge, these features have not been compared in Spiroplasma-infected vs. uninfected embryos.
Our results showing no effect of the male-killing Spiroplasma on maternally-loaded mRNA composition are consistent with what is known about the male-killing mechanism. A Spiroplasma-encoded protein (Spaid) appears to interact with the DCC, which assembles onto the X chromosome of wild-type males but not females, and is accompanied by DNA damage and abnormal apoptosis 15 . This model does not require "priming" of the egg during oogenesis by the symbiont. Nonetheless, the delayed expression of male death in males flies over-expressing wild type Spaid (2nd larval instar), compared to those harbouring male-killing Spiroplasma, requires further investigation and may indicate that other Spiroplasma-encoded factors are relevant to the male-killing phenotype (e.g. RIP) 38 .
The differential expression analyses of the comparisons involving the Hyd1 treatment (i.e., the non-male-killing Spiroplasma strain native to D. hydei), yielded equivocal results, with four pipelines revealing no differentially expressed genes, and the other four pipelines revealing an intersection of five DE genes (ect, Osi6, Osi7, FBgn0038339, and FBgn0037099). All of these five genes had a lower expression level in the Hyd1 treatment compared to the other three treatments (Fig. 4). The time course of expression levels reported in FlyBase for these genes 56 , which is based on Wolbachia-infected flies 57 , indicates that early embryos have lower levels than later embryo stages (see Table 4). Based on this pattern, a plausible hypothesis is that the Hyd1 treatment exhibits a developmental delay. To address this, we compared the expression levels among treatments of six genes that should exhibit high levels at the 0-2 h embryo and substantially lower levels at the 2-4 h embryo (see Table 4). Our rationale was that if the Hyd1 treatment were developmentally delayed, we might detect a trend of higher expression levels at these genes compared to the other treatments. Nonetheless, Hyd1 exhibited no such a pattern; except for one of these genes (i.e., FBgn0035955; not significant; Fig. 7). Therefore, collectively, the patterns of expression of Hyd1-infected embryos are inconsistent with a developmental delay.  www.nature.com/scientificreports www.nature.com/scientificreports/ Available information on the function of the five genes at which the Hyd1 treatment had significantly lower expression provides little insight into the possible causes or consequences. Based on FlyBase, ectodermal (ect; FBgn0000451) is expected to have very low expression (FlyBase RPKM = 1) at the stage we examined. It is expressed later (Stages 13-16 ≈ embryo 14-16 h) in several tissue types: foregut, epidermis, trachea, and hindgut, with possible roles in cuticle development and tubular formation (e.g. tracheal tubes) 58,59 . Osiris 6 (Osi6; FBgn0027527) and Osiris 7 (Osi7; FBgn0037414) are also expected to have very low expression (FlyBase RPKM = 10 and 5, respectively) at the stage examined; with the highest expression at later embryonic stages (i.e., 14-16 h). Both genes are physically close to each other and belong to the Osiris gene cluster (OSI; FBgg0000612). The OSI family of genes shares a highly conserved protein domain is present, and retains high synteny, in all insects examined to date, and appears to have evolved via gene duplications at the base of the Insecta 60,61 . In D. melanogaster, these genes are within the dosage-sensitive Triploid-lethal (Tpl) locus; individuals with one or three copies of Tpl die as late embryos or early first instar larvae 60,62 . Although Osiris genes are generally considered of unknown function, based on expression patterns or knock-down experiments, they have been associated with cuticle formation, metamorphosis, digestion, resistance/tolerance to toxins and viruses, wing formation, and phenotypic plasticity (reviewed in Smith, et al. 63 ). FBgn0038339 (CG6118) contains a BTB/POZ domain, and appears to be involved in regulation of transcription by RNA polymerase II. FBgn0037099 (CG7173) contains a Kazal domain (i.e., a type of serine proteinase inhibitor). Both of these genes are expected to have very low expression (FlyBase RPKM = 1 and 0, respectively) at the stage we examined, with peak expression levels occurring in the 14-16 h embryo.
Assuming that the significant results for the five genes are repeatable, they imply more perturbation of gene expression by Hyd1 than by the other symbiont strains; a finding that could reflect that D. melanogaster is outside  Table 3). The wMel treatment (i.e., Wolbachia-infected) represents the treatment comparable to FlyBase-reported expression values. Nonetheless, the RPKM values reported by FlyBase (and presented in Table 4) are not directly comparable to our normalized counts, for the equivalent treatment.
www.nature.com/scientificreports www.nature.com/scientificreports/ the fundamental niche of Hyd1, as also suggested by its lower vertical transmission efficiency in this host (80% 64 ), and by the difficulty of maintaining this artificial association in the lab (personal observation and Kageyama, et al. 65 ). Reports of additional phenotypic effects of Hyd1 in D. melanogaster have been anecdotal or contradictory. Concerning adult longevity, Kageyama, et al. 65 reported substantially shorter adult life spans of both males and females, whereas Silva, et al. 66 detected little to no adult mortality. The study of Silva, et al. 66 is expected to more closely represent the conditions and specific D. melanogaster (i.e., Canton S) and Hyd1 (i.e., sub-strain TEN104-106) genetic backgrounds of the present transcriptomics study. Concerning fly fecundity, both Kageyama, et al. 65 and Hutchence, et al. 64 report anecdotal evidence of lower female fecundity in Hyd1-infected D. melanogaster than their uninfected counterparts; Silva, et al. 66 did not examine this trait. Consistent with the detrimental effect of Hyd1 on fecundity of D. melanogaster, the microarray-based study of Hutchence, et al. 64 (focused on adults), revealed downregulation of a cluster of genes involved in egg production and fertilization. Nonetheless, in general, Hutchence, et al. 64 found that D. melanogaster adults infected with Hyd1 exhibited less perturbation of gene expression than those infected by MSRO-Br and Spiroplasma poulsonii strain NSRO 64 (a close relative of MSRO that also kills males and naturally occurs in Drosophila nebulosa). The apparently different patterns of gene expression perturbation by Hyd1 between our study and that of Hutchence, et al. 64  The lack of an effect of wMel, a CI-inducing strain, on D. melanogaster early embryo mRNA composition, is not unexpected given what is known about the CI mechanism. Essentially, wMel encodes two contiguous protein coding genes (cifA and cifB) 18 . Each of these proteins is capable of causing CI when expressed in the male germline, whereas only cifA (when expressed in embryos) is able to rescue CI 19 . Our results suggest that wMel does not "prime" the egg for rescue by manipulating maternally-loaded transcripts or their degradation.

Depurination patterns.
Our study capitalized on the presence of (non-target) rRNA reads in RNAseq datasets prepared by poly-A-tail enrichment, to estimate levels of depurination at the RIP target position of the sarcin-ricin loop (SLR) in the 28 SrRNA of eukaryotes. We acknowledge that the levels of depurination estimated from mRNAseq data may be downwardly biased because rRNA depurinated by RIPs is highly prone to hydrolysis of the sugar-phosphate backbone at the lesion site 67 . Furthermore, the freezing to which these samples were subjected between collection and library preparation might have decreased the detectability of depurination 53 . In addition, it is not known whether the polyA-RNA enrichment protocol used prior to library preparation, which is aimed at depleting rRNA in the sample, could bias representation of depurinated vs. intact rRNA. Only one study has compared inferences of ribosome depurination from mRNAseq vs. qPCR assays in the same system. Based on mRNAseq, Hamilton et al. 17 detected ~3.8% depurination (i.e., ~96.2% adenine) in the nematode Howardula infesting adult Drosophila neotestacea harbouring the non-male-killing Spiroplasma sNeo strain. Comparatively, using the qPCR approach, the abundance of depurinated template representing RIP-induced depurination increased ~20-fold, whereas the levels of intact nematode rRNA were reduced ~six-fold in the presence of Spiroplasma 17 . Notwithstanding the potential biases, the consistent finding of no depurination in Spiroplasma-free treatments vs. depurination in Spiroplasma-present treatments across three independent studies, lends credibility to the mRNAseq-based approach we employed for inference of depurination.
The significantly lower proportion of adenines at the RIP target site for the Spiroplasma treatments vs. the Wolbachia treatment and the control is consistent with the action of a Spiroplasma-encoded RIP, and supports recent findings that: (1) MSRO-UG infection causes ribosome depurination in D. melanogaster embryos; (2) that heterologous expression of a Spiroplasma RIP1 and RIP2 genes (separately) also causes ribosome depurination www.nature.com/scientificreports www.nature.com/scientificreports/ in D. melanogaster embryos (measured by a qPCR approach); and (3) that the degree of ribosome depurination ("RIP activity" sensu Garcia-Arraez, et al. 38 ) in the different Spiroplasma treatments or RIP transgene constructs is positively correlated with embryo mortality 38 . Importantly, despite detecting significantly higher levels of depurinated template with the qPCR approach, none of the qPCR assays revealed significantly lower levels of intact template 38 . This suggests that detectable depletion of intact template in the qPCR approach is not a pre-requisite for detecting a phenotypic effect of depurination.
Our results suggest that MSRO substrain Brazil (MSRO-BR) causes comparable levels of depurination to the Uganda substrain, which is consistent with its strong male-killing effect 43,66 , and with the identical content of RIP-encoding genes (unpublished data). Interestingly, the male-killer MSRO generally exhibited a higher signal of depurination than the non-male-killer Hyd1 (MSRO range >1-8%; Hyd1 range < 1-< 3%; Fig. 6). This difference could be the result of differences in titers of MSRO vs. Hyd1 at the embryonic stage. Densities of these strains in D. melanogaster Canton-S background at the adult stage are lower for Hyd1 than for MSRO-Br, but only immediately after adult eclosion 66 ; densities at the embryonic stage have not been compared. Alternatively, it is possible that the RIP genes encoded in the MSRO genome are more actively expressed or secreted, or are more efficient than the putative RIP genes detected in the genome of Hyd1 39 . Levels of fly ribosome depurination in the presence of sNeo, the   Table 4. RPKM counts reported in FlyBase for the five DE genes found in our study and for six additional genes, "Developmental Markers" whose expression is higher in the 0-2 h embryo than in later stages. Gene name or Annotation Symbol in parenthesis. a Note: RPKM values reported in FlyBase are not expected to be directly comparable to the normalized count values reported in the present study. www.nature.com/scientificreports www.nature.com/scientificreports/ non male-killing strain native to D. neotestacea, have not been assayed in embryos, but a qPCR assay of D. neotestacea ovaries indicated a small, albeit non significant, signal of depurination 17 . Similarly, an mRNA-seq experiment of adult D. neotestacea revealed 0.4% depurination (i.e., %A = 99.6%) in the presence of sNeo 17 . Therefore, the three members of the poulsonii clade examined to date (MSRO, Hyd1, and sNeo) encode RIPs capable of depurination of Drosophila ribosomes, but the male-killing strain exhibits the highest levels depurination, a phenomenon that appears to contribute to the male-killing mechanism 38 . The general patterns, however, indicate that Spiroplasma RIPs are particularly efficient at depurinating the ribosomes of natural enemies of Drosophila (i.e., parasitic wasps and nematodes), leading to the hypothesis that RIP-induced depurination plays an important role in defence mechanism 17,37 . The higher depurination levels of fly ribosomes detected in embryos (and ovaries) compared to other fly stages could be due to greater exposure of ribosomes to RIP prior to cellularization, after which Spiroplasma becomes effectively extra-cellular 38 . The relatively high levels of depurination in old adults appears to be the result of higher Spiroplasma densities at that stage 38 . Additional roles have been attributed to RIPs, which could contribute to the male-killing phenotype. For example, several RIPs are reported to cause DNA damage 68,69 . Thus, one or more Spiroplasma-encoded RIP might directly contribute to the DNA damage reported during the process of male-killing 30 .

Conclusions
This study employed a transcriptomics approach to examine whether cytoplasmically-transmitted bacteria, including reproductive manipulators that strongly impact survival of the embryonic stage of Drosophila, influence composition maternally-loaded mRNAs. The results revealed that mRNA composition does not differ significantly among the embryos harbouring the reproductive parasites Spiroplasma MSRO and Wolbachia wMel and those lacking endosymbionts. Only the symbiont Spiroplasma Hyd1, which does not manipulate reproduction, Figure 7. Normalized counts for each treatment and replicate, as estimated by Deseq2, for the six genes regarded as "Developmental markers" because they have high expression in the 0-2 h stage and much lower expression at subsequent stages. FlyBase-reported expression levels (in RPKM) are given in Table 4.