Co‐infection of Wolbachia and Spiroplasma in spider mite Tetranychus truncatus increases male fitness

Abstract Wolbachia and Spiroplasma are intracellular bacteria that are of great interest to entomologists, because of their ability to alter insect host biology in multiple ways. In the spider mite Tetranychus truncatus, co‐infection of Wolbachia and Spiroplasma can induce cytoplasmic incompatibility (CI) and fitness costs; however, little is known about the effect of co‐infection at the genetic level and the molecular mechanisms underlying CI. In this study, we explored the influence of the two symbionts on male mite host fitness and used RNA sequencing to generate the transcriptomes of T. truncatus with four different types of infection. In total, we found symbiont‐infected lines had a higher hatch proportion than the uninfected line, and the development time of the uninfected line was longer than that of the other lines. Co‐infection changed the expression of many genes related to digestion detoxification, reproduction, immunity and oxidation reduction. Our results indicate that co‐infection of Wolbachia and Spiroplasma confers multiple effects on their hosts, and helps illuminate the complex interactions between endosymbionts and arthropods.

Many studies have been performed to elucidate how symbiont infection influences host biology. A deubiquitylating enzyme encoded by Wolbachia was identified to induce CI (Beckmann et al., 2017). In Drosophila melanogaster, prophage WO genes cifA and cifB establish and enhance CI (LePage et al., 2017). In Aedes aegypti and Anopheles gambiae, it was observed that Wolbachia produced striking upregulation of many immune genes and inhibited pathogen development (Kambris et al., 2009;Kambris et al., 2010). In D. melanogaster, Harumoto and Lemaitre (2018) identified a protein from Spiroplasma, SpAID, the expression of which induced male killing. In Drosophila, a ribosome-inactivating protein encoded by Spiroplasma was proven to protect hosts against parasitic wasps and nematodes (Hamilton et al., 2016;Ballinger & Perlman, 2017). However, the biological effects of symbionts on hosts are always determined by host and/or symbiont genotype (Leonardo & Muiru, 2003), and it is necessary to study additional invertebrate systems to reveal the potential mechanisms of reproduction and fitness changes induced by endosymbiont co-infection.
Tetranychus truncatus is a polyphagous agricultural pest that feeds on more than 60 plant species (Bolland et al., 1998) and has become a dominant pest in China in recent years . Field investigation of T. truncatus has shown that co-infection of Wolbachia and Spiroplasma was frequent in nature . Our previous results showed that double infection of Wolbachia and Spiroplasma induced CI in T. truncatus . A general explanation of CI is given by the "modification-rescue model", in which one strain of Wolbachia both modifies sperm during spermatogenesis and "rescues" function to avoid embryonic development disruption (Werren, 1996;Hosokawa et al., 2010). In this model, males play the "modification" role. Although a previous study found that co-infection could confer fitness benefits to their hosts as well as induce CI, little is known about the effect of co-infection at the genetic level and the underlying molecular mechanisms of CI (Werren et al., 2008).
Here, we explored the effects of the two symbionts (Wolbachia and Spiroplasma) on male mite host fitness and compared gene expression profiles among four types of males from different symbiont-infected lines, including doubly infected (Iws), singly Wolbachia-infected (Iw), singly Spiroplasma-infected (Is) and uninfected (U) lines. We found that symbiont-infected lines had a higher hatch proportion than the uninfected line, and the development time of the uninfected line was longer relative to the other lines. Co-infection changed the expression of many genes, including genes related to digestion detoxification, reproduction, immunity and oxidation reduction. Our results indicate that co-infection of Wolbachia and Spiroplasma confers multiple effects on their hosts and helps elucidate the complex interactions between endosymbionts and arthropods.

Mite rearing and sample collection
Tetranychus truncatus were collected from different sites in Shenyang, Liaoning Province, China. Mites were reared on leaves of kidney bean (Phaseolus vulgaris L.) at 25 ±1°C, 60% relative humidity and under 16 : 8 L : D conditions. Three lines, Iws, Iw and Is, were established as previously described . The U line was established by treating the Is line with tetracycline solution (0.1%, w/v) for three generations, and maintaining them in a mass-rearing environment without antibiotic for seven generations before use to avoid potential side effects of antibiotic treatment.
To rule out the effects of host genetic background, females from the Iws, Iw and Is lines were allowed to mate with the males from the U line, and the female offspring continued to mate with males from the U line on new leaf disks. The above-mentioned process was continued for seven generations until identical nuclear backgrounds were obtained. Then, spider mites were cultivated for five generations before experiments.
To ensure experiment accuracy, we checked the symbiont-infected states using polymerase chain reaction (PCR).

Hatchability assay
The effect of symbionts on male mite hatchability was assessed. For each line, a single 1-day-old adult virgin female was placed on a leaf disk and allowed to oviposit for 5 days. Each line used 36 leaf disks. To calculate hatch proportion, the eggs on the leaf disks were checked daily. Data were first tested for normality (Kolmogorov-Smirnov test) and homogeneity of group variances (Levene's test). When the data were normally distributed, they were analyzed with a one-way analysis of variance with post-hoc Tukey Honestly Significant Difference analysis for multiple comparisons. When the data did not follow a normal distribution, Kruskal-Wallis test with Bonferroni correction were used for multiple comparisons. All statistical analyses were done on SPSS 19.0 (IBM, Armonk, NY, USA).

Development assay
To test whether symbiont infection influenced normal male development, we measured development time. For each line, approximately 150 3-day-old adult virgin females were placed on a leaf disk and allowed to oviposit for 3 h (3-day-old adult virgin females were in the active ovipositing state; it increased the possibility of obtaining more eggs in a relative short period). Then, the females were removed and eggs were left. The eggs on the leaf disks were checked at 8:00, 16:00, and 22:00 hours every day to determine whether males underwent eclosion. Observed adult males were picked out and the number of adult males recorded. Development time was considered the time which from an egg was laid to adulthood. The differences of development time were analyzed using the Kruskal-Wallis test with Bonferroni correction. All statistical analyses were done on SPSS 19.0.

Library preparation and RNA sequencing
One-day-old adult virgin males of the four lines were then collected. Samples were stored in liquid nitrogen for subsequent RNA isolation.
Total RNA was extracted with the Trizol protocol according to manufacturer's instructions (Invitrogen, Carlsbad, CA, USA) for a total of 12 independent samples with three biological replicates per treatment. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Libraries were generated using NEBNext R UltraTM RNA Library Prep Kit for Illumina R (NEB, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The complementary DNA (cDNA) library was then sequenced using Illumina Hiseq X Ten (Illumina, San Diego, CA, USA), and paired-end reads were generated.

Quality control
Raw data (raw reads) in fastq format were first processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads that contained adapters, reads that contained poly-N, and lowquality reads from raw data. Simultaneously, Q20, Q30, GC-content and sequence duplication levels of the clean data were calculated. All downstream analyses were conducted with high-quality clean data.

Transcriptome assembly and gene functional annotation
Because there is no reference genome available for T. truncatus, a de novo transcriptome was built with Trinity v2.4.0 (Grabherr et al., 2011), with min kmer cov set to 2 and all other parameters set to default. Gene function was annotated based on the following databases: Nr (National Center for Biotechnology Information [NCBI] nonredundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); EuKaryotic Orthologous Groups (KOG); Swiss-Prot (a manually annotated and reviewed protein sequence database); KO

Quantification of gene expression levels
Gene expression levels were estimated by RSEM (Li & Dewey, 2011) for each sample. Clean data were mapped back onto the assembled transcriptome, and read count for each gene was obtained from the mapping results.

Differential expression analysis
Differential expression analysis of two conditions/groups was performed using the DESeq2 1.6.3 (Love et al., 2014). The P-values were adjusted using the Benjamini and Hochberg (1995) method. Genes with an adjusted P-value <0.05, as determined by DESeq2, were considered differentially expressed.

Enrichment analysis of differentially expressed genes (DEGs)
GO enrichment analysis of the DEGs was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs. KEGG (Kanehisa et al., 2007) is a database resource for understanding highlevel functions and utilities of a biological system, such as a cell, organism or ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). KOBAS (Mao et al., 2005) was used to test the statistical enrichment of DEGs in KEGG pathways.

Quantitative real-time PCR
To estimate the results of RNA sequence (RNA-Seq) analysis, the expression levels of arbitrarily selected genes were measured by quantitative real-time PCR (qPCR). The primer sequences are summarized in Table S1. Total RNA was extracted with the Trizol protocol. cDNA was synthesized according to the manufacturer's instructions (Takara, Kusatsu, Japan). qPCRs were performed on the Applied Bio-systems 7500 Real-Time PCR System with the qPCR SYBR Green Master Mix (Yeasen, China). Each sample was analyzed in triplicate in a 20-μL total reaction volume that contained 0.4 μL (10 μmol/L) of each primer, 7.2 μL ddH 2 O, 10 μL SYBR Green and 2 μL diluted cDNA. The concentration of cDNA was too high for qPCR, so we diluted cDNA 5-fold by adding ddH 2 O in order to get a viable cycle threshold (Ct) value. The relative expression levels were calculated using the 2 − Ct method (Livak & Schmittgen, 2001).

RNA-Seq data processing
Four transcriptome libraries (Iws, Iw, Is and U) were constructed, and each contained 50-54 million reads. A total of 58 015 unigenes were generated from 193 528 396 bp, with a median sequence length of 2408 bp and an N50 of 5468 bp. Over 32 000 unigenes were larger than 2 kb in size ( Table 1). Because of the lack of a T. truncatus reference genome, all unigenes were annotated from the Nr (33827), Nt (11648), KO (21307), SwissProt (34170), PFAM (36542), GO (36868) and KOG (24304) databases (Fig 3). As shown in Figure 4, the vast majority of putative homologous genes were annotated to T. urticae by the Nr database, which was consistent with what we expected.  In total, 661 genes were downregulated and 508 were upregulated in Iws versus U (Fig. 5); 473 genes were downregulated and 546 were upregulated in Iws versus Is (Fig. 6); and 596 genes were downregulated and 514 were upregulated in Iws versus Iw (Fig. 7).

GO analysis
The gene functions of DEGs were annotated using the GO database, and here we list terms which were the most significant enrichment of DEGs in Iws compared with the other three lines. In Iws versus U, DEGs were most enriched in the biological process (BP) categories of amide biosynthetic process and ribosome biogenesis (Fig. 8). In Iws versus Is and Iws versus Iw, fusion of virus membrane with host plasma membrane and membrane fusion involved in viral entry into host cells were the most significant enrichment BP terms. Periplasmic space was the strongest enriched term in the cellular component category in all three comparisons. Additionally, amine dehydrogenase activity and oxidoreductase activity, acting on the CH-NH2 group of donors were the most significant enrichment terms associated with the molecular function category in all three comparisons (Figs. 8,9,10).

KEGG pathway analysis
In Iws versus U, KEGG pathway analysis showed that DEGs were involved in 277 different pathways, mainly including those related to ribosomes and glycolysis/gluconeogenesis (Fig. 11). In Iws versus Is, DEGs were involved in 260 different pathways, mainly including linoleic acid metabolism and arachidonic acid metabolism (Fig. 12). In Iws versus Iw, genes were involved in 263 different pathways, mainly including those related to the phosphatidylinositol signaling system and inositol phosphate metabolism (Fig. 13).

DEGs of interest
Co-infection of Wolbachia and Spiroplasma changed the expression of many more genes compared with other types of symbiont infection ( Table 2).
The polyphagy of T. truncatus could be attributed to powerful detoxification functions, such as having a series of detoxifying enzymes, including cytochrome P450 monooxygenases (CYPs), glutathione S-transferase, ABC transporter, carboxylesterase and multidrug resistance-associated protein. Co-infection of the two symbionts had different degrees of influence on the expression of these genes in a host. For CYP genes, six genes significantly changed in Iws, including four downregulated and two upregulated genes. Four genes that encode glutathione-S-transferase were downregulated (three genes were downregulated compared with Is, and one gene was downregulated compared with U). Additionally, for two genes that encode ABC transporters, one gene was upregulated and one gene was downregulated. Moreover, one carboxylesterase gene was downregulated. Many genes that encode multidrug resistance-associated protein dramatically changed, including five downregulated and six upregulated genes. Several other genes were Fig. 13 The 20 most enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways based on differentially expressed genes between doubly infected (Iws) and singly Wolbachia-infected (Iw) lines. The x-axis shows the rich factor. The y-axis shows the pathway name. The point size represents the number of genes enriched in a particular pathway. Degree of enrichment is more significant for larger rich factor values and smaller q values. also downregulated, including those that encode phosphatase, canalicular multispecific organic anion transporter and P-glycoprotein.
Lipocalins are a family of proteins which transport small hydrophobic molecules. We found that the majority of genes which encode apolipoprotein D were downregulated and only one was upregulated. Several genes associated with larval molt were upregulated. There were several genes related to histone coding and histone modification (Cluster-1760.0, Cluster-6764.36412, Cluster-11188.0, Cluster-6764.8724, Cluster-6764.43656); four of five of these genes were downregulated. We also observed downregulation of three genes related to vitellogenin. Two genes, Cluster-6764.21118 and Cluster-6764.22958, which participate in cuticle biosynthesis, were downregulated. Moreover, it is noteworthy that the genes that encode DNA topoisomerase 2-alpha, which are essential for proper segregation of daughter chromo-somes during mitosis and meiosis, were all substantially upregulated.

qPCR validation
Ten DEGs were arbitrarily picked to confirm sequencing data accuracy (Fig. 14). We examined the expression of these genes in an arbitrary comparison. The differences in expression of all genes except Cluster-6764.13646 were consistent with RNA-Seq results. These results support the reliability of the RNA-Seq data.

Discussion
It is common in nature for multiple endosymbionts to infect the same host (Duron et al., 2008). Co-infection with Wolbachia and Spiroplasma could have various effects    on host reproduction and fitness. Because the effects of symbionts on host biology differ among host and/or symbiont genotypes (Serbus et al., 2008), it is necessary to study additional invertebrate systems to reveal the potential mechanisms of reproduction and fitness changes induced by symbiont co-infection. In addition, a common CI defect is a delay in the breakdown of the male nuclear envelope, which results in paternal chromosome mis-segregation during the first mitosis (Tram & Sullivan, 2002;Tram et al., 2006;Gebiola et al., 2017). Additionally, increased male age can dramatically decrease the strength of CI (Kittayapong et al., 2002;Reynolds & Hoffmann, 2002). Therefore, we collected 1-day-old adult virgin males for the transcriptome analysis in this study.
T. truncatus is a polyphagous herbivore and has numerous genes related to detoxification. In our study, the expression of these specialized genes dramatically changed as a result of symbiont co-infection, which indicates that co-infection of Wolbachia and Spiroplasma affects T. truncatus feeding and digestion. We observed that 34 genes related to digestion and detoxification were differentially expressed in the Iws line, which indicates that the co-infection largely changed host detoxification ability. There were also remarkable changes in many genes which encode canalicular multispecific organic anion transporter and multidrug resistance-associated proteins. They are members of the ABCC subfamily. ABCC proteins have a series of functions, including ion transport, cell-surface receptor activity and translocating many kinds of substrates, such as drugs, cyclic nucleotides, endogenous compounds and their glutathione conjugates and glutathione (Dermauw et al., 2013). The majority of genes that encode cholinesterase were upregulated in Iws. Cholinesterase plays an important role in hydrolyzing both plant secondary metabolites, organophosphates and other pesticides (Ramsey et al., 2010). Many lipocalin genes were differentially expressed. Lipocalins are a family of proteins capable of binding to small hydrophobic molecules (Flower, 1996). In spider mites, they are likely to bind pesticides/allelochemicals, which results in sequestration of these toxic, generally hydrophobic compounds. Furthermore, we found that UDP-glycosyltransferases (UGTs), which catalyze the addition of the glycosyl group from a UTP-sugar to a small hydrophobic molecule, were strongly differentially expressed in Iws. One of the principle functions of insect UGTs is detoxification of ingested plant phenolics, which act as toxins, or insect deterrents (Ahmad & Hopkins, 1993).
In fact, for herbivores, the nutrient content provided by plant material is usually poor or unbalanced (Karban & Baldwin, 2007). The detoxification ability of herbivores enhances utilization of limited plant nutrition. In spider mites, nutritional status is related to host performance (Wermelinger et al., 1991). Our previous study showed that co-infection of Wolbachia and Spiroplasma could confer some fitness benefits to T. truncatus hosts, including shortened development time, and increased fecundity and hatchability, which can all contribute to the expansion of symbionts in host populations . Hence, the detoxification ability of T. truncatus has an indirect impact on the spread of symbionts. We conclude that the strong detoxification ability of the doubly infected line was at least one reason why co-infection of Wolbachia and Spiroplasma is widespread in natural T. truncatus populations Zhang et al., 2016). However, the entirely different symbiont prevalence in the field  probably reflects different levels of detoxification ability of various mite lines. For example, Wolbachia are largely present in the gnathostoma of spider mites , and the function of Wolbachia is related to their localization, which favors the perspective that symbiont infection affects food digestion and detoxification. Further studies are needed to confirm the influence of symbionts on detoxification.
Some genes related to reproduction were differentially expressed as a result of co-infection. Two genes that encode a cuticle protein were downregulated in Iws. In the nematode Brugia malayi, Wolbachia removal resulted in downregulated transcripts involved in cuticle biosynthesis, which may reflect disruption of a normal embryogenic process (Ghedin et al., 2009). Moreover, two genes that encode integrin beta pat-3 were differentially expressed. During Caenorhabditis elegans development, pat-3 is essential for gonad morphogenesis, and pat-3 RNAi caused abnormal gonad formation and sterility (Lee et al., 2001;Lee et al., 2005;Xu et al., 2005). Substantial changes were observed in several genes associated with histone modifications and histone variants. These genes were found to participate in various cell division processes and development. Histone variant H2A.Z plays a critical role in early development, and H2A.Z RNAi results in defective chromosomal segregation (Faast et al., 2001;Rangasamy et al., 2004). Histone methylation is essential for meiosis, because it modulates the expressions of meiotic-related genes (Wang et al., 2017). Additionally, histone acetylation status has a substantial impact on meiosis during both spermatogenesis and oogenesis (Wang et al., 2017).
We observed differential expression of a gene that encodes cyclin-G-associated kinase (GAK). GAK is required for proper mitotic progression because GAK knockdown was shown to cause cell cycle arrest at metaphase (Shimizu et al., 2009). An obvious cytological defect of CI is a delay in chromosome condensation in male pronuclei at prometaphase. Here, we observed differential expression of genes that encode DNA topoisomerase II. DNA topoisomerase II may play a role in chromosome condensation during mitosis (Wang, 2002). Furthermore, it is worth noting that two genes that encode the nuclear receptors ecdysone-inducible protein E75 (E75) and ecdysone-induced protein 78C (E78) were upregulated in Iws compared with Iw. In Tribolium castaneum, E75 regulates male accessory gland development and accessory gland protein production, whereas E78 is associated with sperm production and transfer to females (Xu et al., 2012). In Drosophila, females with low levels of E78 expression experienced a notable decline in egg production and a low hatching rate compared with wildtype individuals (Ables et al., 2015). These genes also regulate larval molting and development time (Bialecki et al., 2002).
In males, development time was found to be is inversely correlated with the strength of CI (Yamada et al., 2007). The spider mites we examined here that were co-infected with Wolbachia and Spiroplasma produced more eggs, and the development time of doubly infected males was relatively short compared with that of males infected with only Wolbachia, and only the Iws line induced CI. Further research is needed to determine the effect of co-infection on the expressions of the above genes in T. truncatus.
Some of the DEGs in the present study that are associated with the immune response were associated with coinfection. These include the genes for the nuclear factor-κB-like transcription factor Relish, a Janus-activated kinase/signal transducers and activators of transcription pathway ortholog of Dm stat92E, Toll-like receptor Tollo and pelle. Moreover, a few genes related to autophagy and apoptosis also dramatically changed (Table S2). We could not confirm that co-infection caused immune system activation, as we did not detect any antimicrobial peptide changes; therefore, further experiments are needed to determine if co-infection activates the immune system. Identifying immune response diversity among different hosts and/or endosymbiont genotypes will also help clarify the interactions between endosymbionts and their hosts.
Furthermore, we observed numerous DEGs related to oxidoreductase activity due to co-infection, such as the genes for dual oxidase, short-chain dehydrogenase and peroxidase, which are important to maintain redox homeostasis.

Conclusions
In summary, we demonstrated that symbiont infection could increase male mite hatchability and accelerate the development from egg to adult. We found that Wolbachia and Spiroplasma co-infection affected transcripts that encode proteins involved in various biological processes in T. truncatus, such as digestion and detoxification, reproduction, the immune system and redox processes.
Our study provides novel gene-based research on the effects of Wolbachia and Spiroplasma co-infection on arthropod hosts, and helps elucidate the complex interactions between symbionts and arthropods. his kind help with the experiment. This study was supported in part by a grant-in-aid for Scientific Research (31672035, 31871976) from the National Natural Science Foundation of China.

Disclosure
The authors declare no conflict of interest.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Primers used for quantitative real-time polymerase chain reaction of differentially expressed genes. Table S2. Genes that were differentially expressed in response to symbiont infection in Tetranychus truncatus males.