Comparative Analysis of Adelphocoris suturalis Jakovlev (Hemiptera: Miridae) Immune Responses to Fungal and Bacterial Pathogens

The wide-spread culture of transgenic Bt cotton resisting the infamous cotton bollworms has reduced the adoption of broad-spectrum insecticides to a large extent. Consequently, the non-targeted insect Adelphocoris suturalis Jakovlev has become a major cotton pest in China. Entomopathogenic microbes show promising results for controlling this pest in the future, but A. suturalis innate immune responses to these pathogens are poorly understood. Here, we used the entomopathogenic fungus Beauveria bassiana and the Gram-negative pathogenic bacteria Enterobactor cloacae to infect A. suturalis nymphs, followed by high throughput RNA-seq to analyze the immune transcriptomes of A. suturalis in response to the two pathogens. A total of 150 immunity-related genes were identified, including pattern recognition receptors, extracellular signal modulators, signal pathways (Toll, IMD, JNK, and JAK/STAT), and response effectors. Further quantitative real-time PCR analysis demonstrated that B. bassiana and E. cloacae were recognized by different receptors (GNBP and PGRP, respectively); activated Toll pathway and IMD pathway respectively; and both induced expression of the effector gene Defensin. However, melanization is suppressed in B. bassiana-infected nymphs. Collectively, this study provides a transcriptomic snapshot of the A. suturalis immune system, and at the genetic level, gains multifaceted insights of the immune response to fungal and Gram-negative bacterial pathogens. Ultimately this work pioneers the study of molecular mechanisms underlying immune interactions between A. suturalis and its pathogens and assists in the development of novel mitigation strategies to control this pest.


INTRODUCTION
Over the past two decades in China, transgenic cottons expressing Bacillus thuringiensis (Bt) toxins have been widely cultured to resist the refractory cotton bollworm Helicoverpa armigera (Wu et al., 2008). This wide-spread control measure substantially reduced the adoption of broad-spectrum insecticides, decreasing selection pressure over Bt-insensitive mirid bug species (Hemiptera: Miridae) . Consequently, one member of mirid bug species, Adelphocoris suturalis Jakovlev, has now become a major cotton pest in China (Li et al., 2010). Biological control methods including transgenic plants, entomopathogenic microbes, or the combination of both, showed promising efficacy in pest's management (Yajuan and Wu, 2010;Luo et al., 2017b). Many studies have focused on A. suturalis ecology and physiology Zhang et al., 2011;Feng et al., 2012;Luo et al., 2017a),reporting its transcriptomes at different developmental stages and pheromone biosynthesisrelated genes (Luo et al., 2014;Tian et al., 2015). Together these reports identified candidate genes for further transgenic plants development. However, the immune responses regarding the interactions between A. suturalis and microbial pathogens, remain poorly understood.
Beauveria bassiana is an entomopathogenic fungus, extensively used as a commercial biopesticide, and which has been proven to competently kill a broad range of pests including Miridae (Boomsma et al., 2014;Mascarin and Jaronski, 2016;Portilla et al., 2019). The infection process starts when a sticky conidium of B. bassiana adheres to the host cuticle, germinates to form a germ tube and an appressorium, and subsequently the cuticle is penetrated by mechanical pressure and hydrolysis of secreted enzymes (Mascarin and Jaronski, 2016). After overcoming the host physical barriers, the hyphae can invade the insect hemocoel, where blastospores proliferate, secreting toxins and depriving the host of essential nutrients, eventually leading to the death of the insects (Boomsma et al., 2014;Valero-Jimenez et al., 2016;Wang and Wang, 2017). Gram-negative bacteria Enterobactor cloacae is a common nosocomial pathogen that can also infect insects and trigger immune responses (Zou et al., 2008;Xiong et al., 2015;Yang et al., 2017). In addition to the potential application in pest management, these two pathogens are also widely used models that comprehensively study insect immune responses (Zou et al., 2010).
In insects, the innate immune system, consisting of physical barriers, humoral, and cellular immune responses, is central to defending against parasitic or pathogenic infections (Lemaitre and Hoffmann, 2007). The physical barrier is mainly attributed to the cuticle (Moret and Moreau, 2012). Successful pathogen invasion activates humoral and cellular innate immune responses. Generally, cellular immunity includes phagocytosis, encapsulation, and nodulation by hemocytes in the hemolymph (Fauvarque and Williams, 2011). Humoral responses result from the production of immune effector molecules and melanization, elicited by the recognition of pathogen associated molecular patterns (PAMPs) by pattern recognition receptors (PRRs) such as peptidoglycan recognition protein (PGRP), Gramnegative binding protein (GNBP), and C-type-lectin (CTL) (Leulier et al., 2003;Gottar et al., 2006). Signal modulation and signal transduction including Toll, IMD, JNK, and JAK/STAT pathways are then subsequently initiated, leading to the generation of effector molecules such as attacin, defensin, and lysozyme (Ferrandon et al., 2007). Where fungal pathogens and Gram-positive bacteria are responsible for the Toll pathway, and Gram-negative bacteria can activate the IMD pathway (Lemaitre et al., 1995(Lemaitre et al., , 1996Michel et al., 2001). The pathogen recognition can also trigger melanization, through activation of prophenoloxidase (PPO), melanin synthesis and pathogen sequestration (Cerenius and Söderhäll, 2004;Cerenius et al., 2008). The insights into immune responses of A. suturalis against different pathogens could improve our understanding of its immune system and guide the development of novel control methods.
In this work, we first investigated the survival rates of A. suturalis injected with B. bassiana and E. cloacae to confirm susceptibility. We then utilized high throughput RNA-seq to analyze the A. suturalis immune transcriptome in response to these two pathogens. The resulting sequence information was mapped to other insect immune-related datasets to identify putative immunity-genes in A. suturalis. Expression patterns of 17 genes from RNA-seq results were further validated by quantitative real-time PCR. This study provided a transcriptomic snapshot of the A. suturalis immune system and gained multifaceted insights of the immune response to fungal and Gram-negative bacterial pathogens at the genetic level. Ultimately this work will pioneer the study of molecular mechanisms underlying immune interactions between A. suturalis and its pathogens and assist in the development of novel mitigation strategies to control this pest.

Insect Rearing and Pathogenic Infection Bioassays
Adelphocoris suturalis nymphs were obtained from Huazhong Agricultural University and reared with fresh mung bean seedlings and pea aphids (Acyrthosiphon pisum) in clean plastic containers at 25 ± 2 • C, 70 ± 5% relative humidity, under the photoperiod of 16 h light/8 h dark (Luo et al., 2017a). E. cloacae (strain no.1.2022) was obtained from China General Microbiological Culture Collection Center and was cultured in LB (Luria-Bertani) medium at 37 • C. B. bassiana isolate 093 was cultured on Potato Dextrose Agar (PDA) medium at 25 • C. As a dispersing agent, 0.05% (v/v) Tween 80 solution was used to collect the fungal conidia. In total, 150 newly molted third-instar nymphs were randomly selected, divided into three parts, and injected with the 0.05% Tween solution (as the control group) or E. cloacae bacterial Tween suspension (about 1.9 × 10 8 cfu mL −1 ), or B. bassiana conidial Tween suspension (about 1 × 10 7 conidia mL −1 ), respectively (Xiong et al., 2015). Each nymph was injected with 20 nL in volume, via a microinjection, using a microinjector (World Precision Instruments, Sarasota, FL, United States). Before the injection, the nymphs were starved for 2 h, and the conjunctives between the metathorax and the first abdominal segment comprised the injection site (Liu et al., 2019). A more detailed process of the microinjection is described in Liu et al. (2019). Survival was counted every 6 h.

cDNA Library Preparation and Illumina Sequencing
To obtain immune transcriptomes with large read coverage, we selected two time points to sample, lethal time of 10% (LT10) and 50% (LT50) (Xiong et al., 2015;Zhang et al., 2017;Xu et al., 2018). Due to the different pathogenicity and lethality dynamics, E. cloacae-infected samples were collected 49 and 84 h post infection. The Tween treatment group was sampled at 49 h after injection, and three live but weak nymphs in each group were pooled as one sample and a whole body of a nymph was sampled here, the total RNA of three biological replications of each treatment were extracted, respectively, so as to extract RNA as much as possible. Total RNA was extracted with an SV total RNA isolation system (Promega, Madison, WI, United States). Then, the mixed RNA of three biological replications of each treatment was used for RNA-seq (Zhang et al., 2015). Five cDNA libraries (Control, Ec_6h, Ec_12h, Bb_49h, and Bb_84h) were constructed using 3 µg of RNA and sequenced on an Illumina Hiseq TM 2000 system (Illumina, United States). FASTQ files of raw reads were produced and sorted by barcodes for subsequent analysis. Raw reads were submitted at the NCBI in the SRA database (PRJNA692314).
All unigenes were annotated against the non-redundant (NR) sequence database and Swiss-prot database using BLASTX (Evalue < 1e −5 ). Gene ontology (GO) annotations for all unigenes were obtained to predict gene functions through the Blast2GO program against the GO database (Ashburner et al., 2000;Conesa et al., 2005). Pathway annotation was performed via the online KAAS tool at the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000).
The RSEM package was utilized to calculate normalized gene expression values, fragments per kilobase per million mapped reads (FPKM) (Trapnell et al., 2010;Li and Dewey, 2011). A differentially expressed gene (DEG) analysis was performed based on the significance level using the R package, DEGseq (p-adjust < 0.001, fold change > 2) (Wang et al., 2010).

Identification of Immunity-Related Genes
Amino acid sequences of known immunity-related genes from other insect species were used as templates for a BLASTX search of A. suturalis unigenes. The immune-related gene database was downloaded from the orthodb database 1 and mapped 1 http://cegg.unige.ch/orthodb7 insect species included Drosophila melanogaster, Anopheles gambiae, Apis mellifera, Bombyx mori, Tribolium castaneum, and so on (Xiong et al., 2015;Zhang et al., 2015). The potential A. suturalis immunity-related genes were identified manually through predictions of similar structure domains against matched genes. The conserved domains were detected using the CD-search tool 2 .

Reverse Transcription-Quantitative PCR Analysis
Reverse transcription-quantitative PCR (RT-qPCR) analysis was adopted to validate immunity-related gene expression patterns from RNA-seq results. RNA samples were prepared using the same method we described above. A total of 2 µg of RNA was reverse transcribed using PrimeScript TM RT Master Mix (Takara, Shiga, Japan) and RT-qPCR was performed on an ABI 7300 Real-Time PCR System (United States) using GoTaq RT-qPCR Master Mix (Promega, Madison, WI, United States). All samples obtained at different timepoints have corresponding Tween controls. Three biological replications were conducted for each treatment. Four technical replications were done for each sample in RT-qPCR. PCR conditions consisted of 94 • C for 5 s, followed by 40 cycles of 59 • C for 20 s and 72 • C for 20 s. The primers used are described in Supplementary Table 1.
The RPS15 gene was used as the reference gene (Luo et al., 2018;Lü et al., 2018). Data were collected, exported to EXCEL, and analyzed by the 2 − Ct method (Schmittgen and Livak, 2008). The mean Ct value of four technical replications in each sample was used for the gene expression analysis, and the relative expression level of each treatment was the average of the three biological replications.

Survival Analysis of Adelphocoris suturalis Nymphs Infected With Beauveria bassiana and Enterobactor cloacae
The survival curves of nymphs infected by E. cloacae and B. bassiana were significantly lower than that of the mock control (nymphs injected with Tween 80, Figures 1A,B, Log-rank test, p < 0.001), demonstrating that A. suturalis is very sensitive to these two pathogens. The LT50 of nymphs infected by E. cloacae was 12 h, and 84 h for nymphs infected by B. bassiana, indicating that E. cloacae has a more efficient lethality to A. suturalis than B. bassiana.

Identification and Functional Classification of DEGs of Adelphocoris suturalis in Response to Beauveria bassiana and Enterobactor cloacae Infection
A total of 29,142 unigenes were annotated based on BLASTX searches, and 5,581 DEGs from the five libraries were obtained according to DEGseq results. Hierarchical clustering of DETs identified five clusters that represented differential expression patterns (Figure 2A and Supplementary Table 2). One-hundredand-twenty-nine genes in cluster 1 were upregulated at 49 h after B. bassiana infection. Transcript levels of 2,776 genes in cluster 2 genes were declined at 6 h but 542 genes in this cluster increased at 12 h post E. cloacae infection. Cluster 3 (553 genes) showed 380 genes increased 6 h after bacterial infection (208 genes increased and 25 genes decreased at 12 h), while 401 genes in cluster 4 showed the gene upregulated 84 h post fungal infection.
According to the Venn diagram analysis of the DEGs, in the B. bassiana-infected nymphs, 1,251 genes were specifically regulated in 49 h (145 up-and 2,795 down-regulated), and 551 genes were exclusively regulated in 84 h (361 up-and 190 down-regulated). In the E. cloacae-infected nymphs, 3,035 genes were specifically regulated 6 h post-infection (228 up-and 2,807 down-regulated), and 596 genes were exclusively regulated 12 h post-infection (295 up-and 301 down-regulated) ( Figure 2B). Based on the GO enrichment analysis, most DEGs in the B. bassiana-and E. cloacae-challenged nymphs were enriched in catalytic activity, binding, metabolic and cellular processes, and in cell parts ( Figure 3A). The KEGG pathway analysis indicated that most DEGs of bacteria-infected nymphs were enriched in ribosome, protein processing in endoplasmic reticulum, carbon metabolism, and purine metabolism (Figure 3B and Supplementary Table 4), and for fungi-infected nymphs, most DEGs were enriched in protein processing in endoplasmic reticulum, RNA transport, and carbon metabolism ( Figure 3B and Supplementary Table 4).

Immunity-Related Genes of Adelphocoris suturalis in Response to Two Pathogens
A total of 150 immunity-related genes were identified ( Figure 4A and Supplementary Table 3) and classified as: pathogen recognition, extracellular signal modulation, intracellular signal transduction (Toll, IMD, JNK, and JAK/STAT pathways), immune response effectors, and others. Thirty-four pathogen recognition receptors were identified, including five scavenger receptors (SR), 12 C-type lectins (CTL), three galectins (GALE), two PGRP, four GNBP, two Thioester-containing proteins (TEP), two FREP genes, two nimrods, and two DRPR (Draper) genes. Eight serine protease inhibitors (Serpin), 20 serine proteases (SP) and three clip domain serine proteases (cSP) involved in signal modulation were identified. One Relish, one Caspar, one Cactus, two Toll, and two spätzle (SPZ) signal transduction molecules were identified in the transcriptome. As for immune response effectors, four prophenoloxidases (PPO) of melanization, three lysozymes (Lys), one lysozyme-like protein (LLP), and one Defensin (DEF) were identified according to the transcriptome analysis. The Venn diagrams of immunity-related DEGs show that fungal and bacterial infections can regulate different immune-related genes ( Figure 4B).

Validation of Immunity-Related Genes Expression Patterns by RT-qPCR
For B. bassiana-infected nymphs at 49 h or 84 h, genes in recognition receptors (GNBP1 and GNBP2), signal modulation (Serpin2, Serpin4, and SPZ), Toll pathway (Toll-1), and effectors (Lys1 and Defensin) were significantly up-regulated, while PGRP2 in recognition receptors and PPO4 in melanization were downregulated. Recognition receptor genes, SP1 and SP7, were downregulated at 49 h post infection and then upregulated at 84 h, while CTL1 showed the opposite pattern ( Figure 4C).

DISCUSSION
In this study, we used the entomopathogenic fungus B. bassiana and the pathogenic bacteria E. cloacae to infect A. suturalis nymphs. Great sensitivity to these two pathogens with high mortality rates was reported (Figure 1). With high-throughput RNA-seq and de novo assembly in the absence of a reference  genome, we matched the pathogen-challenged A. suturalis transcriptomic information with other insects' immune genes, to identify immunity-related genes; and validated expression patterns of some immune genes by RT-qPCR (Supplementary Table 3 and Figure 4C). Altogether this provides a comprehensive insight of the A. suturalis' immune responses to fungal and bacterial pathogens.
Based solely on transcriptomic information from mixed samples, both B. bassiana and E. cloacae sharply evoked the expression of A. suturalis immunity-related genes (Supplementary Table 3). However, B. bassiana and E. cloacae infection triggered differential expression patterns (Figures 4B,C), providing multifaceted information about the A. suturalis immune system. Where for B. bassiana-challenged nymphs, GNBP recognition receptors (GNBP1 and GNBP2) were significantly up-regulated, while PGRP receptors (PGRP1 and PGRP2) were remarkably up-regulated in E. cloacae-infected nymphs (Figure 4C), indicating that different receptors are involved in the recognition of circulating fungal and Gramnegative bacterial pathogens, which is analogous to Drosophila FIGURE 4 | A. suturalis immunity-related genes analysis. (A) Distribution of A. suturalis immunity-related genes in categories of recognition, signaling, regulation, and effectors; (B) Venn diagrams of immunity-related DEGs show that fungal and bacterial infections can regulate different immune-related genes. Overlapping regions show genes that regulated in both samples, while non-overlapping regions represent genes that are specifically regulated in corresponding sample. The directions of transcript level changes are indicated by upward-and downward-pointing arrows; and (C) Quantitative real-time PCR analysis of the A. suturalis immunity-related gene expression, the A. suturalis ribosomal protein 15S gene (RPS15) was used as an internal control. Data are normalized to the expression level in Tween 80-treated nymphs (*p < 0.05). (Lemaitre et al., 1997;Gottar et al., 2002Gottar et al., , 2006. SPZ (Spätzle, SPZ1) and Toll (Toll-1) (Lemaitre et al., 1996) in B. bassianainfected nymphs were up-regulated ( Figure 4C), indicating that this fungus activated the Toll pathway of A. suturalis. Caspar is a suppressor of antibacterial immunity and Relish is a transcription factor of the IMD pathway (Hedengren et al., 1999;Kim et al., 2006). The down-regulation of Caspar and the up-regulation of Relish in E. cloacae-infected nymphs indicated that the IMD pathway was triggered ( Figure 4C). Due to the fact that third-instar nymphs need almost 4 days to develop to fourth-instar nymphs, all the sampled nymphs were still at the third-instar. Theoretically, the nymphs within a same instar, which reared in the same condition with the same treatment, may have similar gene expression patterns. Thus, we chose only one control (Crtl_49h) in the transcriptomic analysis for concision of Figure 2A. Whether the control nymphs at different timepoints have similar gene expression patterns, however, needs to be studied further before being confirmed.
The activation of the Toll and IMD pathways by fungal and gram-negative bacterial infections, respectively, is congruent with previous reports in other insects, such as Drosophila melanogaster (Ferrandon et al., 2007). Defensin, an antimicrobial peptide (AMP) (Richman et al., 1996), was upregulated following infection with both pathogens. Previous work suggested that Defensin is required for the host defense against gram-positive bacteria, but is dispensable against invading Gram-negative bacteria or fungi (Blandin et al., 2002;Ferrandon et al., 2007). We speculate that the B. bassiana and E. cloacae infection may give rise to subsequent opportunistic bacterial infection, but further studies are needed to prove this hypothesis. PPO4 of melanization was down-regulated in B. bassiana-challenged nymphs and up-regulated in E. cloacae-infected nymphs. We speculate that toxins secreted by B. bassiana may restrain the melanization, but again, in-depth studies are required to support such a hypothesis.
Overall, this work provides a comprehensive view of A. suturalis' immune responses to fungal and bacterial infections at the genetic level and will unravel the molecular mechanisms underlying immune interactions between A. suturalis and pathogens. Furthermore, this study will guide the development of efficacious entomopathogenic agents to control this pest.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA; BioProject ID PRJNA692314 and BioSample accession SAMN17319809.