Transcriptomic Immune Response of the Cotton Stainer Dysdercus fasciatus to Experimental Elimination of Vitamin-Supplementing Intestinal Symbionts

The acquisition and vertical transmission of bacterial symbionts plays an important role in insect evolution and ecology. However, the molecular mechanisms underlying the stable maintenance and control of mutualistic bacteria remain poorly understood. The cotton stainer Dysdercus fasciatus harbours the actinobacterial symbionts Coriobacterium glomerans and Gordonibacter sp. in its midgut. The symbionts supplement limiting B vitamins and thereby significantly contribute to the host’s fitness. In this study, we experimentally disrupted the symbionts’ vertical transmission route and performed comparative transcriptomic analyses of genes expressed in the gut of aposymbiotic (symbiont-free) and control individuals to study the host immune response in presence and absence of the mutualists. Annotation of assembled cDNA reads identified a considerable number of genes involved in the innate immune system, including different protein isoforms of several immune effector proteins (specifically i-type lysozyme, defensin, hemiptericin, and pyrrhocoricin), suggesting the possibility for a highly differentiated response towards the complex resident microbial community. Gene expression analyses revealed a constitutive expression of transcripts involved in signal transduction of the main insect immune pathways, but differential expression of certain antimicrobial peptide genes. Specifically, qPCRs confirmed the significant down-regulation of c-type lysozyme and up-regulation of hemiptericin in aposymbiotic individuals. The high expression of c-type lysozyme in symbiont-containing bugs may serve to lyse symbiont cells and thereby harvest B-vitamins that are necessary for subsistence on the deficient diet of Malvales seeds. Our findings suggest a sophisticated host response to perturbation of the symbiotic gut microbiota, indicating that the innate immune system not only plays an important role in combating pathogens, but also serves as a communication interface between host and symbionts.

Introduction which allow the targeted experimental removal of intestinal microorganisms from complex communities are scarce.
Within the Pyrrhocoridae (Hemiptera), characterization of the intestinal microbial community revealed a consistent and conserved microbiota, with the co-occurrence of two actinobacterial taxa belonging to the family Coriobacteriaceae (Coriobacterium glomerans and Gordonibacter sp.) [14,15]. Similar to other hemipteran insects [16,17], firebugs rely on an extracellular posthatch mechanism for the vertical transmission of their actinobacterial symbionts through the deposition of bacteria-containing fecal droplets by adult females over newly laid eggs, with the subsequent probing and uptake of the symbionts by the hatched nymphs [18]. In the European firebug (Pyrrhocoris apterus) and the African cotton stainer (Dysdercus fasciatus) (both Pyrrhocoridae), sterilization of the egg surface results in aposymbiotic individuals, which lack the actinobacterial symbionts [18]. Such individuals show slower growth rates, higher mortality and lower reproductive success [15], which indicates that the actinobacterial symbionts contribute an essential function towards their insect host. A recent study demonstrated that the fitness contribution of the symbionts lies in the nutritional provisioning of B-vitamins to the host, and transcriptomic analyses revealed a tight integration of the symbionts into the host's vitamin metabolism [19].
In this study, we investigated the immune response of D. fasciatus with respect to symbiont perturbation by a comparative transcriptomic analysis of midgut samples from aposymbiotic and symbiotic individuals. We accompany the expression analysis with qPCR validations and phylogenetic analyses of key genes involved in the response to symbiosis such as hemiptericin and defensin. The results of this study provide first insights into the host immune factors involved in the maintenance of this complex and specific extracellular gut microbiota.

Insect rearing and ethics statement
Live specimens of Dysdercus fasciatus were obtained from a laboratory culture at the University of Würzburg, Germany, and reared at the Max Planck Institute for Chemical Ecology, Jena, Germany. The insects were reared in plastic containers (20635622 cm) at a constant temperature of 28˚C and exposed to long light regimes (16 h/8 h light/dark cycles). The bugs were provided ad libitum with previously autoclaved water and dry cotton seeds (Gossypium barbadense).

Experimental elimination of the Coriobacteriaceae symbionts
The generation of aposymbiotic bugs involved the sterilization of egg surfaces following the procedure described previously by Salem et al [15]. Briefly, three day old eggs were submerged in ethanol for 5 minutes, followed by bleach (12% NaOCl) for 45 seconds. Subsequently, residual bleach was removed by washing several times in water. Untreated eggs served as symbiont-containing controls, as it was previously shown that the egg surface sterilization procedure itself does not adversely affect host fitness [15]. To examine the differential pattern of host gene expression in response to symbiont elimination on the natural diet of cotton seeds, five egg clutches of D. fasciatus were harvested, and each clutch was separated into two experimental treatments and reared on cotton seeds: (i) symbiotic (untreated), and (ii) aposymbiotic (egg surface sterilized). This design was chosen to reduce the potential influence of genetic variability among host egg clutches and allowed the use of paired statistical tests for analysis.

RNA extraction and reverse transcription
Three days following adult emergence, a single individual was collected from every experimental treatment replicate, and through dissection, its midgut region (M1-M4) was harvested. Once dissected, the midgut was stabilized in RNAlater solution (Qiagen) and stored at 220˚C. Total RNA isolations were performed using the RNeasy Micro Kit (Qiagen) following the manufacturer's guidelines. Integrity and quality of the RNA samples were determined using the RNA 6000 Nano LabChip kit (Agilent Technologies) on an Agilent 2100 Bioanalyzer (Agilent Technologies) according to the manufacturer's instructions. cDNA was then generated with the QuantiTect Reverse Transcription kit (Qiagen) using the included RT Primer Mix according to the manufacturer's guidelines. To account for possible shortcomings during RNA extraction and reverse transcription and to confirm the success of symbiont elimination after eggs surface sterilization, diagnostic PCR screens targeting the host's 18S rRNA gene and the Coriobacteriaceae symbionts' 16S rRNA genes were conducted using the generated cDNA according to procedures described previously [15].

Transcriptome sequencing and library construction
Sequencing was conducted using RNA extractions from dissected whole midgut regions (M1-M4) of symbiotic and aposymbiotic bugs fed on their natural diet of cottonseeds. RNA from five individuals (one per egg clutch) of each treatment was pooled, resulting in two samples.
Prior to sequencing, the extracted RNA was exposed to rRNA depletion to minimize the amount of ribosomal RNA in the final samples (Analytik Jena, Jena, Germany). Additionally, a poly-A enrichment strategy was used to separately enrich for eukaryotic and prokaryotic RNA, respectively (Analytik Jena, Jena, Germany). Sequencing of the resulting four samples (poly-A enriched and depleted fractions of symbiotic and aposymbiotic bugs, respectively) was performed by a commercial service provider (Fasteris; http://www.fasteris.com) using 5 mg total cDNA ( [20]. The leading and trailing bases of each read were cut off if the quality values were below the default threshold. Additionally, reads were discarded if they were shorter than 30 base pairs after trimming. Following quality checks, the trimmed reads were assembled de novo into contigs using Trinity (r2012-10-05) [21]. The minimal contig length was set to 200 and the kmer length to 25 base pairs. The read libraries of the poly-A depleted and enriched treatments were separately pooled and assembled into backbones, respectively. After the assembly with Trinity, the contigs were clustered by CD-HIT EST (v4.5.7) [22] according to their sequence similarity to remove potential duplicates. Sequences with more than 99% sequence similarity to other contigs were subsequently collapsed. For the assignment of expression values to each constructed transcript in the respective library, the original reads were mapped back to the respective backbone assembly using Bowtie2 (v2.0.3) [23]. The generated output was processed using SAMtools (v0.1.18) [24] to create BAM files and assess the coverage depth as the number of reads mapped to each transcript using the R package Rsamtools (v1.14.1).
The expression information was normalized using the RPKM (reads per kilobase of transcript per million of mapped reads) transformation to obtain estimates of relative expression levels. Homology searches (BLASTx and BLASTn) of unique contig sequences and functional annotation by gene ontology terms (GO; http://www.geneontology.org) were conducted with an E-value cutoff of 10 26 using the BLAST2GO software suite v2.4.1 (http://www.blast2go.de).
Based on the taxonomic annotations of the BLAST search results from the poly-A depleted dataset, we classified different bacterial groups on the order level using MEGAN (v4) [25]. To obtain the relative abundance, we calculated for each bacterial group the overall sum of mapped reads divided by the total number of bacterial reads mapped to the specific group.
In the ploy-A enriched dataset, transcripts with annotations related to the innate immune response were extracted using specific GO terms. Additionally, sequences representing antimicrobial peptides were detected using a bi-directional local BLAST search (E-value cutoff of 10 25 ) with a compiled database consisting of publicly available antimicrobial peptide sequences in Genbank from various insects.

Phylogenetic analyses of selected transcripts
Phylogenetic analyses were performed with the software tool MEGA5 [26] for selected transcripts with the same annotation, as obtained from the NCBI database. The coding regions of the nucleotide sequences were first determined and translated into protein sequences using Augustus (v2.7) [27]. The protein sequences were subsequently aligned with Muscle [28]. Maximum likelihood trees were then constructed with 500 bootstrap replicates.

Validation of host gene expression by quantitative PCR
Quantitative PCRs (qPCRs) for the candidate host immune effector genes were conducted across the individual RNA samples from aposymbiotic and symbiotic bugs used for sequencing (five replicate individuals per treatment) to confirm the transcriptome sequencing results. Primers were designed based on the candidate gene sequences of hemiptericin, defensin, c-and i-type lysozyme (S1 Table) from the transcriptome and checked for specificity in vitro using capillary sequencing of amplified PCR products as described previously [15]. For pyrrhocoricin, no sufficiently specific and efficient qPCR assay could be designed, based on the short transcript sequence available, so this gene was not included in the qPCR analysis. Conditions for qPCR were optimized using a VWR Gradient Thermocycler (VWR, Radnor, PA, USA) at various annealing temperatures (60-68˚C).
The qPCR reactions were performed using a RotorGene-Q cycler (Qiagen, Hilden, Germany). The final reaction volume of 25 ml included the following components: 1 ml of cDNA template, 2.5 ml of each primer (10 mM), 6.5 ml of autoclaved distilled H 2 O, and 12.5 ml of SYBR Green Mix (Qiagen, Hilden, Germany). Standard curves for absolute quantification in the qPCR (10-fold dilution series from 1 ng/ml to 10 26 ng/ml) were generated using purified PCR products for all primer pairs after measuring the PCR product concentrations using a NanoDrop1000 spectrophotometer (Peqlab, Erlangen, Germany). For qPCR, the following cycling parameters were used: 95˚C for 10 min, followed by 45 cycles of 68˚C for 30 s, 72˚C for 20 s, and 95˚C for 15 s. Subsequently, a melting curve analysis was conducted by increasing the temperature from 60˚C to 95˚C within 20 min. Six replicates of one of the standard concentrations were used, for each primer pair and concentration, for the configuration and calibration of the standard curve. The resulting averages were then utilized to correct for possible errors in the DNA concentration measurements. Based on the standard curve, absolute transcript copy numbers were calculated according to [29].

Statistical analysis
The absolute expression determined by qPCR was normalized with the expression of a housekeeping gene (60S ribosomal protein L13a). The normalized expression levels in the aposymbiotic and symbiotic treatments were then compared using the Wilcoxon signed rank test to assess levels of significance (P,0.05).

Success of symbiont elimination
Sterilization of egg surfaces resulted in aposymbiotic firebugs that were free of both Coriobacteriaceae symbionts, as confirmed by C. glomerans-and Gordonibacter-specific diagnostic PCRs. Conversely, symbiotic bugs tested positive for both bacterial species (data not shown).

Bacterial transcripts in the midgut of D. fasciatus
In total, RNA sequencing of the poly-A depleted treatments yielded 38,207,274 reads with a combined length of 3.8 Gbp. After trimming, backbone assembly and subsequent clustering of similar sequences resulted in 89,807 contigs with a N50 length of 955 base pairs. Despite the poly-A separation strategy, only a small subset of bacterial transcripts (1,651 contigs) were recovered and annotated in the poly-A depleted fraction (Fig. 1), which, in turn, failed to provide a coherent overview of the microbiota's metabolic and cellular processes in D. fasciatus' midgut. However, the transcriptional signal was still indicative of the taxonomic diversity of metabolically active bacterial species within the midgut of symbiotic and aposymbiotic firebugs.
Overall, the expression patterns between aposymbiotic and symbiotic bugs were highly similar with minor differences in Lactobacillales, Rickettsiales and Rhizobiales. The strongest difference was found in relation to transcripts belonging to the Coriobacteriaceae symbionts (C. glomerans and Gordonibacter sp.). As expected and in contrast to the high expression patterns observed in symbiotic bugs, no Coriobacteriaceae transcripts were retrieved from the aposymbiotic treatment.
Immune system-related transcripts in the midgut transcriptome of D. fasciatus In total, RNA sequencing of the poly-A enriched treatments yielded 36,456,544 reads with a combined length of 3.6 Gbp. After trimming, backbone assembly and subsequent clustering of similar sequences resulted in 55,222 contigs with a N50 length of 1,411 base pairs. Of those contigs, 20,174 received an annotation via BLAST homology. The majority of the annotations were based on hits from Riptortus pedestris and Acyrthosiphon pisum, the bean bug and pea aphid, respectively. From the BLAST results, 3,362 contigs were assigned to functional categories according to gene ontology (GO), most of which were determined to be involved in cellular and metabolic processes and biological regulation.
Here we focus on gene candidates involved in immune-system processes of the insect host (Fig. 2), recovered from the poly-A enriched treatment. Four genes putatively responsible for recognition, and 20 genes involved in signal transduction in the major insect immune pathways Toll, Imd, JAK/STAT, and phenoloxidase were recovered (Table 1, Fig. 2). With nine annotated transcripts, the Toll pathway constituted the most comprehensive set, whereas the Imd and JAK/STAT pathways contained significantly fewer annotated transcripts (both four transcripts). For the phenoloxidase pathway, three transcripts were assigned as pro-phenoloxidase (Dfas-51289), serine protease easter-like (Dfas-14165), and mature phenoloxidase (Dfas-45148).
The expression of transcripts presumably involved in recognition and signaling generally showed relatively low expression levels based on the normalized RPKM values ( Table 1). The low expression pattern of these transcripts nonetheless pointed towards an up-regulation (Fold .2) in aposymbiotic individuals (Table 1), with the exception of the transcripts annotated as protein toll (Dfas-01522), nf-kappa b inhibitor (Dfas-48512) and peptidoglycan recognition protein s2 (Dfas-25404), which were up-regulated (Fold .2) in symbiotic bugs. Most of the transcripts potentially involved in signaling showed no differential expression at all.

Antimicrobial peptide transcripts
Among the immune effector genes, several AMP types with potential isoforms were detected (Table 2)  33105), two as i-type lysozyme (Dfas-45802, Dfas-30083), and one as c-type lysozyme (Dfas-30397). The immune effector genes generally showed much higher expression levels than genes putatively involved in immune signaling and microbe recognition (Table 1 and 2). Comparison of AMP expression between symbiotic and aposybiotic bugs revealed that most of the detected transcripts were differentially expressed (Fold .2), with some genes being up-and others downregulated in response to symbiont elimination ( Table 2). Among the transcripts up-regulated in aposymbionts were i-type lysozyme, hemiptericin, and defensin, whereas pyrrhocoricin and c-type lysozyme were down-regulated in symbiontdeprived bugs. Within each AMP, the expression of the isoforms was relatively consistent, with only minor variations in the RPKM values.

Expression validation of antimicrobial peptides with qPCR
To validate the observed normalized expression patterns of c-and i-type lysozyme, defensin and hemiptericin (Table 2), quantitative PCRs were performed (Fig. 3). The same expression patterns were observed for hemiptericin and c-type lysozyme compared to the RNAseq expression analysis. The expression of

Phylogenetic analysis of transcripts encoding antimicrobial peptides
To study the micro-diversity of sequences annotated as defensin and hemiptericin, phylogenetic analyses were performed (Fig. 4). For both genes, sequence alignments revealed various amino acid substitutions across the different isoforms (S1 and S2 Figures). Defensin-related transcripts showed two distinct isoforms, one of which clustered with a published sequence for P. apterus (Fig. 4A). For hemiptericin, phylogenetic analysis revealed three distinct isoforms (Fig. 4B), one of which clustered together with sequences of P. apterus obtained from the GenBank database. The transcript annotated as c-type lysozyme ( Table 2) was highly similar to the R. pedestris lysozyme sequence (S3 Figure).

Discussion
In this study, we conducted a comparative transcriptomic analysis of midgut samples from symbiont-containing and aposymbiotic individuals of the cotton stainer, D. fasciatus. We identified a number of differentially expressed transcripts involved in the immune response of the host following symbiont elimination (Fig. 2). These results provide first insights into the molecular interactions between cotton stainers and their extracellular intestinal symbionts and complement the knowledge of immune responses in perturbed gut ecosystems of insects [30]. In the following paragraphs, we will discuss the putative biological roles of the immune related genes identified in this study and the implications of their differential expression upon manipulation of the microbial midgut community in D. fasciatus.

Effect of symbiont elimination on the gut microbiota
The few collected transcripts of bacterial taxa (Fig. 1) were insufficient to provide a comprehensive view of the microbial gene expression. However, they did provide significant insights into the taxonomic composition of metabolically active bacterial symbionts in the midgut of D. fasciatus. Our analysis revealed a microbial community that is largely consistent with earlier results characterizing the gut microbiota of firebugs, with high abundances of Actinobacteria, Firmicutes, and Proteobacteria [14,15]. Most importantly, bacteria of the family Coriobacteriaceae were abundant in symbiotic bugs, but completely absent from aposymbiotic individuals. These results augment previous findings [15] demonstrating that the egg surface sterilization is specific towards ridding firebugs from their association with C. glomerans and Gordonibacter sp.
Concerning the distribution of other microbial taxa, only minor quantitative differences between aposymbiotic and symbiotic bugs were detected (Fig. 1). This indicates that even though a major group of mutualistic bacteria were removed from the intestinal ecosystem, the host was still able to maintain a specific composition of the remaining microbiota, by providing a selective environment and/or by regulation via the immune system. The minor differences in gene expression of transcripts belonging to Lactobacillales, Rickettsiales and Rhizobiales might be in this respect indicative for a compensatory mechanism from the host side.

Constitutive expression of immune signaling pathways
We identified a subset of the known pattern recognition proteins (Table 1) involved in the sensing of foreign organisms [8]. Among those was one short type peptidoglycan recognition protein (PGRP s2-like). In Drosophila, short class peptidoglycan recognition proteins are responsible for the detection of Grampositive bacteria and subsequent induction of immune pathways [31]. Thus, this protein may be involved in the recognition of the mutualistic Coriobacteriaceae in D. fasciatus, which is consistent with the up-regulation in symbiotic as compared to aposymbiotic bugs (Table 1).
Among the immune-related transcripts detected in this study (Table 1, Fig. 2), all four major insect immune pathways [8] were recovered in D. fasciatus. The detected Toll pathway provided the most comprehensive view, with annotated transcripts spanning the signal transduction cascade from precursor protein spaetzle to the transcription factors dorsal and cactus [32]. In contrast, the detected genes involved in the Imd pathway provide only a fragmented view of this signaling cascade. The scarcity of annotated Imd-pathway genes is likely due to the absence of many Imd genes in the pea aphid Acyrthosiphon pisum [33], which is the closest organism to the Pyrrhocoridae with a sequenced genome available, and thus, with reliable homology information. Across the four major pathways (Toll, Imd, JAK/STAT, phenoloxidase), gene expression was generally low. Furthermore, no distinct regulation patterns could be detected, with only few genes being differentially expressed (.2 fold) between symbiotic and aposymbiotic individuals. These observations are in line with other studies, demonstrating the constitutive expression of immune signaling pathway genes [34].

Expression patterns of immune effectors depict differentiated biological roles
In contrast to the low and constitutive expression of pathway-internal genes, the transcripts of the end products of signal transduction, the immune effector genes, show generally higher expression levels that differ between symbiotic and aposymbiotic bugs ( Table 2). Considering that antimicrobial peptides (AMPs) are the front line of defense against microbial pathogens [9], this indicates a highly specific as well as regulated immune response towards symbiont perturbation.
The expression patterns across the detected antimicrobial peptides differed considerably (Table 2). In general, two groups could be differentiated: i) Transcripts down-regulated in aposymbiotic bugs comprised c-type lysozyme and pyrrhocoricin, whereas ii) defensin, hemiptericin and i-type lysozyme were upregulated upon symbiont elimination. Validation of expression patterns via qPCR (Fig. 3) confirmed the differential expression for c-type lysozyme and hemiptericin, but not for defensin and i-type lysozyme across the studied individuals. Therefore, c-type lysozyme and hemiptericin were considered strong candidates for an important role in responding towards perturbation of the microbial midgut community. Additionally, pyrrhocoricin represents an interesting candidate gene with potential down-regulation in aposymbiotic bugs, but the lack of a specific qPCR assay prevented us from gaining more detailed insights on the consistency of the expression patterns across replicate individuals.
The high expression of c-type lysozyme (Fig. 3) and pyrrhocoricin (Table 2) in symbiotic bugs might serve to regulate the association with the pivotal Coriobacteriaceae symbionts. This idea is supported by the primary activity of both lysozymes [35] and pyrrhocoricin [36] against Gram-positive bacteria. In weevils of the genus Sitophilus, the expression of a single antimicrobial peptide gene (coleoptericin A [colA]) was found to be important in controlling population dynamics of mutualistic bacteria, which overreplicated in the host, when colA was knocked down with RNA interference [6]. The same might be true in Pyrrhocoridae species, which may require c-type lysozyme to control the abundance of the intestinal symbionts in order to maintain the host-symbiont equilibrium. This equilibrium may be particularly important for gut ecosystems, since at a higher microbial growth rate, the host might compete with its symbionts for important nutrients rather than engage in a mutualistic relationship.
Furthermore, considering the nutritional relevance of the symbiotic Coriobacteriaceae [19] up-regulation of c-type lysozyme might cause the release of nutrients via lysis of the bacterial cells [19]. In Drosophila melanogaster, lysozyme expression plays an important role in the digestion of intestinal bacteria [37], indicating a similar pattern as observed in this study (Fig. 3). This harvesting of nutrients through lysis of symbiont cells was also proposed in R. pedestris as a mechanism to release nutrients and control symbiont populations [38].
In contrast to c-type lysozyme and pyrrhocoricin, hemiptericin was significantly up-regulated in aposymbiotic individuals, and defensins as well as itype lysozyme showed tendencies towards higher expression in aposymbiotic bugs. As in other insects, hemiptericin is known to be important in P. apterus for the response to foreign microorganisms, particularly Gram-negative bacteria [36]. The up-regulation upon symbiont elimination in D. fasciatus might therefore constitute a compensatory mechanism to control the remaining gut microbiota, which might otherwise overreplicate in the ecological niche freed by the absence of the Coriobacteriaceae.

Isoforms of antimicrobial peptides indicate a global response to the gut microbiota
Interestingly, several isoforms of hemiptericin and defensin were detected (Fig. 4), which may allow for a fine-tuned and dynamic interaction with the intestinal community, since the different isoforms have divergent sequences (S1 and S2 Figures) and may thus act against different microbial taxa. This is particularly relevant in cotton stainers, as in addition to the Coriobacteriaceae symbionts, they consistently harbor several bacterial taxa in the Firmicutes and Proteobacteria (Fig. 1). These taxa could be differentially regulated in their abundance by individual AMPs, with the defensins primarily targeting Gram-positive and hemiptericins Gram-negative bacteria [36]. The presence of AMP isoforms has also been described for other insects [13,39], but their role and adaptive significance in immunity remains unknown.

Conclusions and perspectives
Using a comparative transcriptomics approach, our study revealed a differentiated immune response to the elimination of a specific group of extracellular gut symbionts. This response includes both up-and down-regulation of certain antimicrobial peptides of the host, respectively, which may serve to maintain a stable equilibrium of the complex microbiota in the host's gut and could additionally be involved in harvesting nutrients through lysis of symbiont cells. While the transcriptomic approach provides a global view into gene expression patterns in the cotton stainer's gut, targeted knock-down of individual immune effector genes by RNA interference [40] is necessary to provide functional insights into the interactions between the microbial symbionts and the host immune system.