First transcriptome of the Neotropical pest Euschistus heros (Hemiptera: Pentatomidae) with dissection of its siRNA machinery

Over the past few years, the use of RNA interference (RNAi) for insect pest management has attracted considerable interest in academia and industry as a pest-specific and environment-friendly strategy for pest control. For the success of this technique, the presence of core RNAi genes and a functional silencing machinery is essential. Therefore, the aim of this study was to test whether the Neotropical brown stinkbug Euschistus heros has the main RNAi core genes and whether the supply of dsRNA could generate an efficient gene silencing response. To do this, total mRNA of all developmental stages was sequenced on an Illumina platform, followed by a de novo assembly, gene annotation and RNAi-related gene identification. Once RNAi-related genes were identified, nuclease activities in hemolymph were investigated through an ex vivo assay. To test the functionality of the siRNA machinery, E. heros adults were microinjected with ~28 ng per mg of insect of a dsRNA targeting the V-ATPase-A gene. Mortality, relative transcript levels of V-ATPase-A, and the expression of the genes involved in the siRNA machinery, Dicer-2 (DCR-2) and Argonaute 2 (AGO-2), were analyzed. Transcriptome sequencing generated more than 126 million sequenced reads, and these were annotated in approximately 80,000 contigs. The search of RNAi-related genes resulted in 47 genes involved in the three major RNAi pathways, with the absence of sid-like homologous. Although ex vivo incubation of dsRNA in E. heros hemolymph showed rapid degradation, there was 35% mortality at 4 days after treatment and a significant reduction in V-ATPase-A gene expression. These results indicated that although sid-like genes are lacking, the dsRNA uptake mechanism was very efficient. Also, 2-fold and 4-fold overexpression of DCR-2 and AGO-2, respectively, after dsRNA supply indicated the activation of the siRNA machinery. Consequently, E. heros has proven to be sensitive to RNAi upon injection of dsRNA into its hemocoel. We believe that this finding together with a publically available transcriptome and the validation of a responsive RNAi machinery provide a starting point for future field applications against one of the most important soybean pests in South America.

Identification of RNAi-related genes. The result of the E. heros transcriptome search for RNAi-related genes revealed the presence of 47 genes associated with dsRNA uptake, RNAi core machinery, auxiliary RISC factors, nucleases, antiviral RNAi, and intracellular transport. Some RNAi-related proteins presented variants, with the presence or absence of conserved domains. Overall, the sequences of H. halys showed the highest similarity to sequences from E. heros. dsRNA uptake. The protein sequences involved in dsRNA uptake were searched in the E. heros transcriptome, and a total of six proteins related to this process was found, although there was an absence of sid-like genes ( Table 1, Supplementary data S1). Scavenger protein was found with a CD36 domain region, Ubiquitin-protein transferase (FBX011) with an F-box conserved domain and three beta-helices, and Epsin 2 with an Epsin N-terminal homology (ENTH) domain. The Clathrin heavy chain (Chc) protein and Gap Junction Protein with an Innexin conserved domain were also found in the E. heros transcriptome.
Core RNAi machinery. Proteins related to the miRNA, siRNA and piRNA pathways were identified in the E. heros transcriptome ( Table 2, Supplementary data S2).
The DCR-1 protein was found in E. heros with the conserved PAZ (Piwi, Argonaute and Zwille) domain, two RNaseIII domains and a Double-stranded RNA-binding domain (DSRBD), with an absence of the helicase domains. DCR-2 was also found in E. heros with two isoforms as following: 1 and 2 with 646,601 and 0.618 transcripts per million (TPM), respectively. The DCR-2 isoform 1 contained all the conserved domains: one helicase domain, one PAZ domain, two RNaseIII domains, and a DSRBD, while DCR-2 isoform 2 was found with two RNaseIII domains and a Ribonuclease III C terminal domain (RIBOc). Dicer 3 protein was not found in the E. heros transcriptome. Drosha protein was found with two RNaseIII domains and a RIBOc, but with the absence   www.nature.com/scientificreports www.nature.com/scientificreports/ identified: AGO-1, AGO-2, AGO-3, Aubergine (Aub) and Piwi ( Table 2, Supplementary data S2). Four variants of the AGO-1 protein were found: isoforms 1, 3, 4 and 5, presenting 0.585, 0.103, 0.4 and 118,437 TPM, respectively. All AGO-1 isoforms were found with the PAZ and PIWI conserved domains. For the AGO-2 protein, two isoforms were found, isoform 1 and 2, with 146,222 and 0.14 TPM, respectively. AGO-2 isoform 1 was found with PAZ and PIWI conserved domains, while AGO-2 isoform 2 had no PAZ domains. AGO-3, Aub and Piwi proteins presented the PAZ and PIWI conserved domains. Zucchini (Zuc), with a nuclease conserved domain, was also found in the E. heros transcriptome.
Phylogenetic analyses revealed distinct groups for DCR and AGO superfamily proteins ( Supplementary  Fig. S1, S2). The protein DCR-1 from E. heros was grouped in a clade with the DCR-1 proteins from Nezara viridula (Hemiptera: Pentatomidae) and H. halys, and the same results were found for E. heros DCR-2 ( Supplementary  Fig. S1). Also, E. heros DCR-1 was grouped in a distinct clade compared to E. heros DCR-2, but it showed a common ancestor. The phylogenetic analysis of the AGO superfamily resulted in two main clades, one contained the AGO subfamily proteins, AGO-1 and AGO-2, while the other had the PIWI subfamily proteins, AGO-3, Aub and Piwi ( Supplementary Fig. S2). E. heros AGO-1 was clustered with AGO-1 from N. viridula and N. lugens, while E. heros AGO-2 was clustered in a second group together with other proteins of this family. E. heros AGO-3 was clustered in a distinct group as well as Aub and Piwi proteins.
Auxiliary RISC factor. The E. heros transcriptome was searched for RNAi auxiliary factors (Table 3, Supplementary data S3). The research resulted in 17 intracellular factors associated with the RISC. The Tudor-SN (TSN) protein sequence, with a Tudor-conserved domain, and the Translin and Translin-associated factor-X (TRAX), conserved subunits of the component 3 promoter of the RISC (C3PO), were identified in E. heros. The Armitage (Armi), spindle-E (Spn-E), Maelstrom, Gawky, Staufen (STAU) and CLIP-associating protein  Nucleases. Exoribonuclease 1 (Eri-1) and DNA/RNA non-specific endonuclease (dsRNase) proteins were found in the E. heros transcriptome ( Ex vivo dsRnA hemolymph degradation. The dsRNA stability in the hemolymph was assessed at 0, 1, 10, 30, 60 and 120 min of incubation. After 10 min of incubation, the dsRNA-V-ATP-A was partially degraded, as the gel showed a smear below the band, clearly demonstrating dsRNA degradation ( Supplementary Fig. S5). At 30 and 60 min of incubation we observed increased degradation, with all dsRNA degraded after 120 min incubation.  www.nature.com/scientificreports www.nature.com/scientificreports/ Mortality of E. heros by dsRNA microinjection. Mortality was assessed at 24, 48, 72 and 96 h post-microinjection of dsRNA-V-ATP-A ( Fig. 3). At 24 h, there was 7% mortality and this increased to 19% at 48 h, 28% at 72 h, and at 96 h 35% of the treated insects were killed. Alongside the mortality in dsRNA-V-ATP-A treated E. heros, reduced mobility was observed compared to the insects microinjected with dsRNA-GFP, which were very active. These mobility effects lasted until 72 h., with a recovery in the mobility at 96 h post microinjection.

Gene expression of V-ATPase-A, and DCR-2 and AGO-2 in E. heros. The V-ATPase-A transcripts
level gradually decreased following dsRNA treatment over time (14% to 74% from 24 h to 48 h, respectively) ( Fig. 4). At 72 h and 96 h, there was an increase in the relative transcript levels, with ~40% reduction in gene expression, but despite this, these values were still significantly lower than the control (dsRNA-GFP-microinjected) insects (p-values <0.001 and 0.014, respectively) (Fig. 4).
The involvement of the siRNA machinery in the gene silencing mechanism was assessed through a qRT-PCR analysis of DCR-2 and AGO-2 genes (Fig. 5). The relative transcript levels of DCR-2 were significantly higher in the insects microinjected with dsRNA-V-ATP-A compared to the controls (not exposed to dsRNA), with the highest DCR-2 expression level observed at 72 h post microinjection and with an increase of ~2.0-fold ( Fig. 5A). At 96 h, relative transcript levels of DCR-2 dropped to ~1.5-fold, still higher than the controls. The expression pattern of AGO-2 behaved similarly as we saw for DCR-2 ( Fig. 5B). At 48 and 72 h post-microinjection, the relative transcript levels of AGO-2 were higher with almost a 4.0-fold increase compared to the control samples.

Discussion
The Neotropical stinkbug E. heros is one of the most important soybean pests in Brazil, Argentina, and Paraguay, and the current lack of genetic information is among the factors limiting the prospects of RNAi as an alternative control approach.
The RNAi pathway works primarily through dsRNA uptake, intracellular dsRNA transport, dsRNA processing to sRNA, RISC complex formation and binding, and digestion/repression of the target mRNA 48   www.nature.com/scientificreports www.nature.com/scientificreports/ our currently reported E. heros transcriptome database, most of the genes involved in these processes above and related to RNAi pathways, are also present in the E. heros transcriptome (Tables 1-4). However, it is important to note that although these genes are involved in the RNAi process in other organisms, it does not mean that they play the same role in the RNAi mechanism in E. heros, and the real involvement of these genes needs to be further confirmed in future functional assays.
To achieve gene silencing through RNAi, dsRNA is taken up by the tissue/cell. In eukaryote organisms, this process occurs through sid-like transmembrane proteins 25,49 or endocytosis-mediated uptake 24,25 . Before Sid-like homologous proteins have been found in Coleoptera, Hymenoptera, Lepidoptera, and Hemiptera, but not in Diptera 49 . Also in the E. heros transcriptome, sid-like homolog genes were not found. In Drosophila (Diptera: Drosophilidae), with the lack of sid-like homolog genes, the dsRNA uptake occurs via endocytosis-involving scavenger receptors 23,24 . Indeed previous work demonstrated that the Scavenger protein is involved in endocytic dsRNA uptake in insects 23,24,50 and other organisms, such as mites 19,51,52 . The Chc protein, which is related to an alternative mechanism for endocytic dsRNA uptake in insects [23][24][25][26] , was found in the E. heros transcriptome. Consequently, with the absence of sid-like genes in the E. heros transcriptome, we believe that the Chc protein may be involved in cellular uptake in E. heros; however the involvement of this protein in dsRNA uptake needs to be proven in future functional assays. In addition, future experiments need to investigate the importance of endocytosis in E. heros.
Core RNAi machinery genes were also searched for in the E. heros transcriptome with focus on the miRNA, siRNA and piRNA pathways (Table 2), and most of these were present with the absence of a RNA-dependent RNA polymerase (RdRP) gene. The lack of RdRP was generally expected because, so far, it has been reported only in ticks, plants and in the nematode Caenorhabditis elegans (Rhabditida: Rhabditidae) 53 . The main core domains of Dicer are well known due to their involvement in dsRNA cleavage into small RNA molecules (sRNAs), including miR-NAs and siRNAs. In the current work, the DCR-1 protein, which is related to the miRNA pathway, contains a PAZ domain, two RNaseIII domains, and a DSRBD, however no conserved domains for helicase were identified. For DCR-2, two isoforms with distinct structures and abundances were identified. The DCR-2 isoform 1 was the most abundant and showed a helicase domain, a PAZ domain, two RNaseIII domains, and a DSRBD 54,55 . The PAZ domain holds a binding pocket for the 3′ overhang of dsRNA substrate and a phosphate-binding pocket that recognizes the phosphorylated 5′ end of small RNAs 56,57 . The two RNaseIII domains are the catalytic core components of Dicer and responsible for the cleavage of the dsRNA substrate 58 . The function of the helicase domain remains unclear, but so  www.nature.com/scientificreports www.nature.com/scientificreports/ far it is known that this domain is required to process siRNA but not miRNA 56 . In flies, the loss in the functionality of DEAD/Helicase domain is related to a particular function in the miRNA-based gene regulation 56 . We hypothesize here that the loss of the helicase domain in the DCR-1 protein in E. heros may be a functional adaptation, related to the miRNA pathway, but this needs to be further investigated. The canonical conserved domains of DCR-2 isoform were not identified in E. heros. Similarly, DCR-2 isoforms with the lack of conserved domains were also identified in mammals 59,60 as well as in Arabidopsis thaliana 61 . Due to the lack of important functional domains it is expected that these DCR-2 variants may not be involved in the siRNA pathway, however the function of these isoforms in insects still remains unclear. To our knowledge, this is the first report of DCR-2 variants in insects. It would be interesting in the future to investigate the role of DCR variants in cellular processes.
In other insects such as Cylas puncticollis (Coleoptera: Curculionidae), N. lugens, D. v. virgifera, L. decemlineata, Drosophila and Tribolium, DCR-1 and DCR-2 are also present 31,48,62,63 . In Drosophila the involvement of DCR-1 and DCR-2 is well established in the miRNA and siRNA pathways 62 . In the piRNA pathway, there is no evidence of a dsRNA precursor and the need of DCR endonucleases [64][65][66] . Drosha protein was identified with two RNaseIII domains plus a RIBOc 54 and with some similar features to Dicer, although it processes miRNA precursors in the nucleus 17 . The dsRNA-binding proteins Pasha, Loquacious and R2D2, which mediate dsRNA binding to the RISC complex, are among the other proteins from the DCR superfamily identified in E. heros. These proteins are cofactors required to interact with the RNaseIII genes Drosha, DCR-1, and DCR-2, respectively 31,62 (Table 2).
Five members belonging to the Argonaute superfamily were identified in the E. heros transcriptome as follows: AGO-1 and AGO-2 which belong to Argonaute subfamily, and AGO-3, Aubergine (Aub) and Piwi, which belong to the PIWI subfamily 67,68 . AGO-1 is an essential protein related to the miRNA pathway, and AGO-2 is related to the siRNA pathway 68,69 . More recently, two new functions have been attributed to AGO-1 and AGO-2 in early Drosophila melanogaster embryos: the generation of polarity within cells and tissues by modulating an important cell-cell signaling pathway 70 . These proteins are characterized by the presence of PAZ and PIWI domains, which guide sRNA recognition and binding, supporting endonucleolytic cleavage 71 . The PAZ domain forms a pocket for siRNAs binding and, specifically, the characteristic two nucleotides (nt) 3′ overhangs, trimmed by Dicer proteins, while the PIWI domain shares structural similarities with ribonucleases and degrades the corresponding RNAs [72][73][74] . The lack of a PAZ functional domain in the AGO-2 isoform 2 raises the hypothesis that this isoform may be related in another biological process, as mentioned above for the DCR-2 isoforms. In the shrimp Marsupenaeus japonicus (Decapoda: Penaeidae), three AGO-1 isoforms have been identified, and interestingly, two isoforms were more expressed in the lymphoid organ, suggesting a role in immunity 75 . The presence of multiple isoforms of AGO-1 and AGO-2 may indicate a role of AGO in many biological processes, including cell proliferation/ differentiation, immune defense, among others 75 . AGO-3, Aub, and Piwi are proteins related to the piRNAs pathway 65,68 and they were also found in E. heros (Table 2). Zucchini (Zuc), responsible for piRNA maturation 76 and related to the germline RNAi processes 77 , was also identified in the E. heros transcriptome.
The identification of both DCR-1 and DCR-2 was confirmed through a phylogenetic analysis using sequences from other insects and revealed distinct groups inside DCR proteins (Supplementary Fig. S1). The E. heros protein DCR-1 was grouped in a clade with DCR-1 proteins from N. viridula and H. halys, showing a common ancestor. The same results were found for E. heros DCR-2. Also, E. heros DCR-1 was grouped in a distinct clade compared to E. heros DCR-2, but with a common ancestor. The phylogenetic analysis for the AGO superfamily resulted in two main clades; one containing the AGO subfamily proteins AGO-1 and AGO-2, and another with the PIWI subfamily proteins AGO-3, Aub and Piwi (Supplementary Fig. S2). E. heros AGO-1 clustered with AGO-1 from N. viridula and N. lugens with the same ancestor. E. heros AGO-2 was assembled in a second group together with other proteins of this family. These two clusters showed a common ancestor. The AGO-3 was clustered in a distinct group, as well as Aub and Piwi proteins. E. heros AGO-3 was grouped with the AGO-3 proteins from H. halys and other insects. Thereby, the phylogenetic analyses were useful to confirm the identification of the core RNAi genes present in the E. heros transcriptome. www.nature.com/scientificreports www.nature.com/scientificreports/ AGO protein is the core component of the RISC, and guided by the siRNA it promotes mRNA cleavage 72,73 . Next to AGO, other important genes related to RISC were identified in the E. heros transcriptome (Table 3). Tudor-SN (TSN) protein is known to interact with Argonaute proteins in the silkworm Bombyx mori (Lepidoptera: Bombycidae) 78 , while Translin and TRAX, that are conserved subunits of the C3PO, are involved in RISC activation, supporting RNAi activity 79 . The Armi, Spn-E, Maelstrom and Hen-1 nucleases are involved in piRNA biogenesis 31 . Maelstrom mutations in Drosophila ovaries resulted in a depletion of Dicer and AGO-2 proteins, the latter two being related to the RNAi pathways 80 . Elp-1, that is also present in E. heros, is a component of the polymerase II elongator complex, and although the absence of this protein in Drosophila S2 cell lines did not affect the miRNA pathway, it can cause an inhibition of the siRNA pathway 81 . The Vasa intronic gene (VIG), that encodes a putative RNA-binding protein through association with RISC 82 , and related to the production of piRNAs 83 , was also identified in the E. heros transcriptome. The Gawky protein, a cytoplasmic mRNA component necessary in early embryonic development 84 , Staufen (STAU), a DSRBP, and CLIP-associating protein (Clp-1), that is responsible for the phosphorylation of the 5′ end of siRNAs 85 and related to the splicing process of transfer RNAs 86 , were all also found in the E. heros transcriptome of this study. The PRP16 protein plays a role in the pre-mRNA processing 87 , while Belle has a function in the endo-siRNA pathway 88 . The proteins GLD-1 and ACO-1, known to inhibit translation of mRNA into protein 81 , were also identified in E. heros.
Nucleases (RNases together with other RNA enzymes) function in DNA/RNA digestion in the midgut 89 and offer an additional defense and regulatory control layer. The activity of nucleases in dsRNA degradation (dsR-Nases) is well known, taken an important role in RNAi efficiency across insect groups such as Hemiptera 39,41,90 , Lepidoptera 91,92 , and Diptera 38 . Four nucleases were identified in the E. heros transcriptome: Eri-1, Nibbler, SDN1, and dsRNase (Table 4; Supplementary Fig. S3). The Eri-1 nuclease is suggested to play a role in the intracellular siRNA and miRNA pathways 93 . In C. elegans, Eri-1 forms a complex with Dicer, generating specific classes of siRNAs, while in mouse, Eri-1 negatively regulates the global abundance of miRNA 93 . Nibbler, an exonuclease known to be involved in shaping the 3′end of the miRNAs, and its depletion leading to developmental defects in Drosophila 94 , was found with conserved domains in E. heros. Another intracellular nuclease found in E. heros was SDN1. In Arabidopsis, this protein is involved in the degradation of mature miRNA, and the knockdown resulted in developmental defects 94 . However, the involvement of the Eri-1, Nibbler, and SDN1 in RNAi efficiency in insects remains unclear. In E. heros we also identified a dsRNase gene with seven isoforms and with a conserved Endonuclease_NS domain associated with the degradation of foreign dsRNA molecules 38 . In B. mori three dsR-Nases isoforms were identified and expressed in different tissues, such as epidermis, fat body, and gut; these dsRNases are related to the innate immune response against invasive nucleic acids 92 . The presence of a dsRNase nuclease with seven isoforms may indicate that E. heros has a strong nuclease activity, so this may result in a lower potential to suppress the expression of target genes and so in turn a lower RNAi response.
In the current work, we identified some genes related to antiviral RNAi as follows: Ars2, a gene related to RISC regulation, ninaC, a gene associated with vesicle transport, and a seven transmembrane-domain glycosyltransferase, egh 47 (Table 4; Supplementary Fig. S4). These genes are known to be involved in antiviral defense in Drosophila 47,95 . The CG4572 gene was also identified in E. heros; it is a carboxypeptidase with unknown function, but related to RNAi in D. melanogaster 47 . Three genes involved in intracellular transport were also identified. Two vacuolar H+ adenosine triphosphatases (V-ATPases) genes were identified in the E. heros transcriptome: V-ATPase subunit A (vha68) and V-ATPase subunit C (vha16). These genes are located at different functional V-ATPase domains, the peripheral domain (V1) and the integral domain (V0) 96 , respectively, and they are related to dsRNA release by the endocytic vesicles 23 . The Small Rab GTPases and vha68 are essential signaling components linked to the extracellular part with the cytoplasm in L. decemlineata 25,48 .
The presence of some genes in E. heros suggests that it has an active and functional RNAi machinery. However, the lack of sid-like gene and the presence of nuclease raise the concern about the RNAi efficiency. So, we first checked the stability of a dsRNA molecule in the hemolymph of adults in which it was rapidly degraded. After 10 min, the dsRNA-V-ATP-A was partially degraded, with increasing dsRNA degradation over time up to 120 min ( Supplementary Fig. S5B). In a similar experiment with the pea aphid A. pisum, dsRNA was completely degraded after 3 h incubation and this was associated with the lack of RNAi responses in this species 41 . In E. heros, dsRNA was completely degraded after 2 h of incubation with watery saliva 97 . Indeed, high nuclease activity in the hemolymph and saliva of E. heros may reduce RNAi efficiency and so some form of dsRNA protection may be needed for future field applications.
To confirm the effectivity of the E. heros RNAi machinery, a dsRNA targeting the V-ATPase-A gene was microinjected into adults. Previously, targeting the V-ATPase-A gene led to mortality in Pectinophora gossypiella (Lepidoptera: Gelechiidae) 98 , E. heros nymphs 97 , A. pisum 99 , H. halys nymphs 42 , among others. The main V-ATPase function is the pumping of protons across the membrane 100,101 , generating an energy gradient. E. heros adults were microinjected with ~28 ng of dsRNA-V-ATPase-A per mg of insect fresh weight, and the mortality was evaluated at 24, 48, 72 and 96 h after microinjection. At 24 h post-microinjection, there was 7% mortality and this increased to 35% at 96 h (Fig. 3). The same dsRNA concentration previously demonstrated to cause up to 50% mortality in E. heros 2 nd instar nymphs 7 days post-microinjection 97 . Based on these results, we believe that this species is sensitive to RNAi when we compare to other insects also considered sensitive to RNAi. Fishilevich and collaborators 102 used the same dsRNA concentration through microinjection in E. heros adults targeting chromatin remodeling genes and this significantly reduced fecundity and egg viability. Coleoptera insects are considered to be more sensitive to RNAi, presenting a robust RNAi mechanism, while Lepidoptera and Hemiptera appear to be more recalcitrant 103 . Second-instar larvae of the African sweet potato weevil C. puncticollis were microinjected with 200 ng/mg of body weight targeting different genes, and mortality reached up to ~50% after six days 104 ; this concentration is ~9 times higher than that one used in E. heros. One of the main reasons associated with the lack of RNAi response in the C. puncticollis weevil was the high nuclease activity 104 www.nature.com/scientificreports www.nature.com/scientificreports/ 26% at 96 h post-microinjection in 3 rd -instar larvae 98 . One strategy to increase RNAi efficiency is an adequate formulation of the dsRNA molecules. In E. heros nymphs, liposome-encapsulated dsRNA targeting V-ATPase-A led to 45% mortality after 14 days as compared to 30% with naked dsRNA 97 . Similar results were found for dsRNA α-tubulin and lipoplexes in the German cockroach, Blattella germanica (Blattodea: Blattellidae) 106 . Therefore, the formulation of dsRNA may provide an affordable non-transformative easy-to-use strategy to deliver gene silencing for pest control in the field. However, for successful pest control, it is very important to know the dsRNA concentration, expressed as per mg of insect body weight, to permit a rationalized pest control strategy based on dsRNA concentration and the delivery approach.
Alongside the mortality, other effects were also observed. The treated insects exhibited reduced mobility in contrast to the insects microinjected with dsRNA-GFP which were very active. This effect lasted until 72 h post-microinjection. Retardation in larval development was reported in P. gossypiella 98 and Helicoverpa armigera (Lepidoptera: Noctuidae) 107 treated with dsRNA targeting the V-ATPase.
To confirm that the observed mortality in E. heros injected with dsRNA-V-ATPase-A is a true phenotype of gene silencing, a qRT-PCR assay was performed. Indeed we confirmed the V-ATPase-A gene silencing with a reduction of 74% in the relative level of transcripts. At 72 and 96 h post-microinjection, there is an increase in gene expression but still 40% lower compared to the insects treated with GFP (Fig. 4). An increase in DCR-2 and AGO-2 gene expression was observed with the highest gene expression observed at 48 and 72 h, with a respective increase of ~2.0 and ~4.0-fold, so confirming the activation of the siRNA machinery upon exogenous dsRNA delivery (Fig. 5A,B). This data has shown the activity of the siRNA machinery in E. heros through the supply of dsRNA. As expected, due to the high nuclease activity, the RNAi effects were temporary, and at 72 and 96 h, there was a recovery in the relative transcript levels from the target genes. Also, at 96 h post-microinjection, there was a reduction in the expression of the siRNA-related genes. In Manduca sexta (Lepidoptera: Sphingidae), an upregulation of DCR-2 and AGO-2 expression in response to injection with dsRNA with also only transient effects in the gene upregulation was reported 108 . As discussed above, DCR-2 and AGO-2 are core components of the siRNA pathway, and the overexpression of DCR-2 and AGO-2 after dsRNA microinjection confirmed the upregulation of the siRNA machinery.
Over the past year, scientists have made enormous progress towards the use of RNAi as a pest control strategy taking advantage of genetic sequences available in public databases, and used this information to understand the RNAi mechanism in insects. To our knowledge, this is the first study of E. heros transcriptome, including the identification of RNAi-related genes and dissecting the siRNA pathway. The analyses of the E. heros transcriptome have identified the main components of the three RNAi pathways with the surprising lack of sid-like genes. Identification of the core RNAi genes, efficient mortality rates and activation of the siRNA machinery, these data provide a novel and important dataset on RNAi machinery and its efficiency, underpinning future strategies to enhance RNAi in E. heros and potentially other piercing-sucking insects as models or species important in agriculture.

Material and Methods
Brown stink bug insects. The colony of E. heros was originally started with insects collected in Pelotas, Brazil (27°48′1.7352″ S; 52°54′3.834″W) in 2013, and kept for about 73 generations under laboratory conditions before experiments. New insects collected in soybean fields in Rondinha, Brazil (27°48′1.7352″ S; 52°54′3.834″ W) were introduced in the laboratory colony in 2015. All stages were maintained in plastic cages under laboratory conditions with a photoperiod of 14:10 (Light: Dark), temperature of 25 ± 1 °C and 75 ± 10% relative humidity. Green beans, peanut and water were supplied ad libitum and replaced twice in a week. Eggs were collected twice a week to obtain the insects necessary for microinjection and colony maintenance 109 . Insects were collected every day and insects of four days old were used in the microinjection assays. cDNA libraries, Illumina sequencing and de novo assembly. Eggs, all nymphal stages and adults of E. heros were used for total RNA extraction using the Trizol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer instructions. The RNA pool was prepared with an equally RNA amount from all stages, and the cDNA library preparation and Illumina sequencing were conducted at the Laboratory of Functional Genomics Applied to Agriculture and Agri-Energy, at the University of São Paulo, Brazil. The TruSeq RNA Sample Prep kit (Illumina) protocol was used to construct the cDNA library, following manufacturer instructions. A high-throughput Illumina sequencing platform (HiSeq. 2000) was used for the final library sequencing, in one lane of a 100 bp paired-end run.
The raw reads originating from the Illumina sequencing were check for quality using the FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). After that, reads were trimmed using Trimmomatic 110 , and only high-quality reads, showing a Phred score superior to 30, were used for the de novo assembly to generate a set of contigs using Trinity software (http://trinityrnaseq.sourceforge.net) 111 . De Bruijn graph algorithm and a k-mer length of 25 were used as parameters.
Homology search and gene ontology annotation. The generated contigs were analyzed using the UniProt-TrEMBL database 112 via Diamond algorithm 113 , with an E-value<10 −5 as a cut-off parameter. The contigs with insect hits were submitted to a second homology search using QuickGO to identify gene ontology (GO) terms. For this annotation, a similarity search was performed against the UniProt database using Diamond, with an E-value<10 −5 as a cut-off parameter. The QuickGo from EBI (https://www.ebi.ac.uk/QuickGO/annotations) was used to calculate the GO terms, using the generated gene identifiers as inputs.
RNAi-related genes. We searched for the genes related to RNAi efficacy and these included genes on dsRNA uptake (Table 1), RNAi core machinery (Table 2), auxiliary factors (Table 3), nucleases, antiviral RNAi (2020) 10:4856 | https://doi.org/10.1038/s41598-020-60078-3 www.nature.com/scientificreports www.nature.com/scientificreports/ and intracellular transport (Table 4) 31,33,48 . Homologous sequences for these proteins were searched in UniProt or Protein database from NCBI, and we used as a query to search the E. heros transcriptome using the tBLASTn tool from NCBI. Generated contigs with a bitscore>150 and E-value <1e-5 were further used to confirm the identity. To detect the open reading frames (ORFs) in the contigs sequence, the ORF Finder from NCBI was used (https://www.ncbi.nlm.nih.gov/orffinder/), and the protein domains predicted by the NCBI Conserved Domains Database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). Protein Basic Local Alignment Tool (Protein BLAST) was used for protein homology search against insect non-redundant protein database at NCBI.
To provide additional confirmation on identity and function prediction of the core RNAi proteins, nucleases and antiviral RNAi, members of these groups of proteins were subject to a phylogenetic analysis using the neighbor-joining (MEGA 7.0.26) algorithm with 1,000 bootstrap replicates. A total of 39 Argonaute superfamily protein sequences, 30 endoribonuclease III protein sequences, 25 nuclease protein sequences, and 27 antiviral RNAi sequences were aligned using the MUSCLE program from MEGA 7.0.26 software. ORF Finder from NCBI was used to predict the proteins. dsRNA synthesis and purification. Specific primers were used to amplify the fragments of the target genes ( Table 5). The cDNA was synthesized using the SuperScript First-Strand Synthesis System Kit (Invitrogen) following the manufacturer's instructions. The T7 primer sequence (TAATACGACTCACTATAGGGAGA) was placed in the front of the forward and reverse primers. These primers were used for dsRNA synthesis with cDNA as a template. The PCR reaction was performed with 2 µl of cDNA template, 2 µl of a 10 µM solution of each primer (Integrated DNA Technologies, Coralville, IA, USA), 0.125 µl of Taq DNA polymerase, 2.5 µl of Buffer 10×X, 0.5 µl of 10 µM dNTPs, 0.75 µl of MgCl 2 (Invitrogen) and 15.125 µl of nuclease-free water (GE Healthcare, Little Chalfont, UK) in a total volume of 25 µl. The PCR conditions used were 5 min at 94 °C for initial denaturation, followed by 30 s at 94 °C, 45 s at 59.5 °C, 55 s at 72 °C for 30 cycles and final extension for 10 min at 72 °C. The amplified products were purified using a PCR purification kit (Qiagen, Valencia, CA, USA) and analyzed on 1% agarose gels. The PCR product was quantified using a Nanovue spectrophotometer (GE Healthcare) and then samples were stored at −20 °C.
The V-ATPase-A dsRNA was synthesized using the MEGAscript T7 RNAi kit (Ambion, Austin, TX, USA) following the manufacturer's instructions. The control group consisted of a dsRNA of the green fluorescent protein (dsRNA-GFP) synthesized from a DNA plasmid (pIG1783f) and cloned in Escherichia coli (DH5α). Plasmid DNA was extracted and sequenced to confirm the identity of PCR products. The identity of the sequence was confirmed by Sanger sequencing. The dsRNA was analyzed for integrity on 1% agarose gels, its concentration quantified in a Nanovue spectrophotometer (GE Healthcare) and then stored at −20 °C.
Ex vivo dsRnA hemolymph degradation assay. Insects were anesthetized with CO 2 during ~30 s and then taped with the abdomen upwards on a glass plate. Legs and rostrum were cut, and hemolymph collected by a needle, prepared with glass capillary tubes, coupled to an insulin syringe (8 ×X 0.30 mm) and placed in chilled 1.5 ml tubes containing phenylthiourea (PTU) to prevent melanisation. After that, 30 µl of dsRNAs-V-ATPase-A solution at 200 ng/µl was incubated in 3 µl of RNase-free water or 3 µl of hemolymph at 25 °C. Aliquots of 5 µl were collected after 0, 1, 10, 30, 60, and 120 min, and the same volume of EDTA (10 mM) was added to the solution to stop the enzymatic reaction 41 . The integrity of the samples was analyzed on 1% agarose gel.

Adult microinjection.
To silence the V-ATPase-A gene in E. heros, dsRNA-V-ATP-A with 623 bp was microinjected in adults (~60 mg) at the concentration of ~28 ng per mg of body weight (0.50 µl of a 3350 ng/µl dsRNA solution) 102 . The control group consisted of insects microinjected with a 560 bp dsRNA molecule targeting GFP 31,104 . The dsRNA-V-ATP-A was designed to have a length similar to the one used in previous RNAi assays in the hemipterans D. citri 43 and N. lugens 114 .
To perform the microinjection, insects were anesthetized with CO 2 and immobilized in a glass plate with double-sided tape (3 M, São Paulo, Brazil). The microinjection was performed using an insulin syringe (8 ×x 0.30 mm) with a needle (30 g) (Solidor) coupled to a micro-applicator (Burkard, Rickmansworth, UK). In total, 62 adults were injected per treatment, of which 12 individuals were used for qRT-PCR and 50 individuals for mortality assay, at 24, 48, 72 and 96 h post-microinjection. Alongside mortality analysis, visual observations were carried out to analyze other effects related to the dsRNA in the insects. After microinjection, the insects were placed in plastic cages containing green beans, peanut and water ad libitum, and kept at 25 °C, photoperiod of 14:10 (L:D) and 75 ± 10% RH, as with the colony maintenance. The insect mortality was normalized against the control (dsRNA-GFP).
Quantitative reverse-transcription polymerase chain reaction (qRT-PCR). Total RNA was extracted from whole insect body at 24, 48, 72 and 96 h after microinjection, and each time point had three biological samples containing one insect. RNA extraction was performed using RNAzol RT (MCR, Cincinnati, OH, USA), following the manufacturer's instructions. The samples were quantified using a Nanovue spectrophotometer (GE Healthcare), verified in a 1% agarose gel electrophoresis, and kept at −80 °C. First-strand cDNA synthesis proceeded as described in the dsRNA synthesis and purification section.
The qRT-PCR was performed on a Roche LightCycler 480 (LC480) (Roche Diagnostics, Basel, Switzerland) real-time PCR platform. To validate the primers used in the analysis (Table 5), a melting curve analysis with temperatures from 60 to 95 °C and a standard curve based on a serial dilution of cDNA were used to determine the primer annealing efficiency and specificity. The reaction included 6 µl of EvaGreen 2X qPCR MasterMix (ABM, Milton, ON, Canada), 1.25 µl of 10 µM forward primer (Integrated DNA Technologies), 1.25 µl of 10 µM reverse primer (Integrated DNA Technologies), 2.5 µl of nuclease-free water and 2 µl of cDNA, in a total volume of 13 µl. The amplification conditions were 3 min at 95 °C followed by 45 Table 5. Primers used in qRT-PCR and dsRNA synthesis.