RNA interference in the Asian Longhorned Beetle:Identification of Key RNAi Genes and Reference Genes for RT-qPCR

Asian Longhorned Beetle (ALB) Anoplophora glabripennis is a serious invasive forest pest in several countries including the United States, Canada, and Europe. RNA interference (RNAi) technology is being developed as a novel method for pest management. Here, we identified the ALB core RNAi genes including those coding for Dicer, Argonaute, and double-stranded RNA-binding proteins (dsRBP) as well as for proteins involved in dsRNA transport and the systemic RNAi. We also compared expression of six potential reference genes that could be used to normalize gene expression and selected gapdh and rpl32 as the most reliable genes among different tissues and stages of ALB. Injection of double-stranded RNA (dsRNA) targeting gene coding for inhibitor of apoptosis (IAP) into larvae and adults resulted in a significant knockdown of this gene and caused the death of 90% of the larvae and 100% of adults. No mortality of both larvae and adults injected with dsRNA targeting gene coding for green fluorescence protein (GFP, as a negative control) was observed. These data suggest that functional RNAi machinery exists in ALB and a potential RNAi-based method could be developed for controlling this insect.

insects. Several other factors such as the digestion of dsRNA by dsRNase enzymes present in the gut lumen and hemolymph, the dsRNA transport into and within the cells and its systemic spread and even the expression levels of RNAi genes have been shown to influence RNAi efficiency among insect tested [12][13][14][15][16][17] .
To determine the knockdown efficiency of dsRNA, reverse transcription quantitative PCR (RT-qPCR) is often used. Although RT-qPCR is a robust method with a combination of high specificity, accuracy, and sensitive detection, this method is influenced by several factors, such as the stability of reference genes, quantity and purity of RNA used, and primer efficiency 18,19 . A critical step that may compensate most of the variability at both RNA and PCR level is the selection of a reference gene constitutively expressed across the samples and treatments under study 20 . Recent studies have suggested that there is no universal reference gene that is appropriate to normalize gene expression in different organisms and conditions 21,22 . To facilitate the validation of reference genes, four models based on different statistical algorithms, Genorm 20 , NormFinder 23 , BestKeeper 24 and delta-Ct 25 were combined in a free web tool, RefFinder 26 . RefFinder provides an overall final ranking based on the calculation of the geometric mean of each program to estimate the stability of the candidate reference genes 26 .
The current studies were designed to identify and validate the RNAi machinery of A. glabripennis, as well as to identify reliable reference genes for normalization in gene expression studies. In silico analysis was performed to identify genes involved in RNAi pathways, based on their homology to sequences of RNAi genes identified in other insects. This was then followed by validation of a reliable reference gene and the profiling of iap (inhibitor of apoptosis) gene expression. RNAi bioassays were performed with both ALB larvae and adults using dsRNA targeting iap gene. Injection of dsRNA targeting iap (dsIAP) into larvae and adults caused knockdown of target gene and death of both larvae and adults. This is the first study to provide evidence for RNAi response and validation of reference genes for expression analyses for ALB. These data could help future studies on gene expression as well as the development of RNAi-based methods against this and other invasive wood-boring insects.

Results and Discussion
Identification of ALB RNAi genes. A total of 33 RNAi genes were identified in the genome of ALB (Table 1). The genes coding for proteins involved in siRNA (Dicer 2, R2D2, Ago 2), miRNA (Dicer 1, Drosha, Loqs, Pasha, Ago 1), and piRNA (Aubergine, Ago 3) pathways were identified. The Dicer proteins contain two RNAse III, helicase, Dicer, Paz, and dsRNA binding domains. Argonaute proteins contain PAZ and PiWi domains. Drosha contains two RNase III domains. Loquacious, R2D2, and Pasha contain three, two and one dsRNA binding domains, respectively. A phylogenetic tree was constructed by comparing the sequences genes identified in ALB with their homologs in other organisms (Fig. 1).
The presence of genes coding for proteins involved in the transport of dsRNA into and within the cells as well as to spread the silencing signal throughout the organism is very important for successful RNAi. Discovered in C. elegans, systemic RNA interference deficient (SID) proteins are essential and sufficient to mediate uptake and systemic spread of RNAi signal 27 . While SID-1 homologous genes have been identified in some insects, SID-2 and SID-5 28 , and the tyrosine kinase SID-3 29 have not been reported in insects. In ALB, we identified two Sid1-like proteins (Sil) homologous to SilA and SilC in Tribolium castaneum. An additional Sil (SilB) protein was identified in T. castaneum and Bombix mori 11 . However, we did not find SilB homolog in ALB. Whether or not SilA and SilC proteins are involved in dsRNA transport in ALB requires further investigation. We also identified homologs of T. castaneum genes coding for proteins known to function in vesicle-mediated transport (Arf72A, AP50, and Rab7), intracellular transport (IdICP, Light, and NinaC), as well as Rsd-3, Lqf (orthologous to Epn-1/Lqf), and scavenger receptors in ALB. Some of these proteins could act as dsRNA receptors and function in dsRNA transport and systemic RNAi mechanism 11 . Other proteins identified to function in the endocytosis pathway in Leptinotarsa decemlineata, such as clathrin heavy chain and vacuolar H + ATPase the 16 kDa subunit, subunit c (Vha16) 30 , were also identified in this study. These data showed the presence of RNAi machinery in ALB. In addition, we identified three dsRNases homologous to dsRNase genes in other insects 31 (Table 1). In Lygus lineolaris (Hemiptera) dsRNases identified in saliva and completely digest long dsRNAs resulting in no-silence effects by orally administered dsRNA 32 . Further investigation of AgladsRNases including determining the expression levels of these two genes in specific tissues of this insect may uncover the role of these enzymes in the RNAi efficiency.
Validation of reference genes. For normalization of gene expression in RT-qPCR, a moderately expressed reference gene is preferred because extremely high or low expression of a housekeeping gene could introduce variability in data analysis, so a standard Ct value range was analyzed for all three experiments (Fig. 2). The expression levels of all candidate genes were measured by the Ct value, which is the number of PCR cycles needed to reach a specific threshold level of detection and is inversely correlated with the quantity of initial RNA template. Succinate dehydrogenase flavoprotein subunit A (sdf) and Ubiquitin (ubq) showed lower levels of expression (Ct value 24). Tubulin (tubulin), ribosomal protein L32 (rpl32), glyceraldehyde 3-phosphate dehydrogenase (gapdh) and elongation factor-1 alpha (ef1a) showed moderate levels of expression (Ct values between 16 to 20 cycles).
Four different programs were used for analysis of reference gene expression (geNorm, NormFinder, BestKeeper, and delta-Ct method) to estimate the stability of six candidate reference genes among different tissues and development stages using a web tool that provides a reference gene stability ranking. A final ranking based on the calculation of the geometric mean of the four algorithms was generated by RefFinder, where the smaller the geometric mean, the greater the stability of reference gene expression. The first experiment compared gene expression among different larval tissues. The geNorm program was used to calculate the stability of the reference genes based on an M value. The lower the M value, the more stable is the expression of the reference gene, and values of M that surpass the cutoff value of 1.5 are not considered stable across treatments. According to this algorithm, all candidate genes had M ≥ 1.5. The gapdh, ubq and ef1a genes showed the lower M value of 1.03, and consequently the most stable genes. The NormFinder program analyzes both intra and inter-group variations, and lower output scores indicate the reduced variation of the reference gene expression. The gapdh, ubq and ef1a are the most stable reference genes with M value of 0.79, 0.89, and 0.97. The BestKeeper algorithm calculates standard deviation (SD), with lower values considered more stable, and values that surpass the cutoff value (SD < 1) are considered to be unstable across all treatments. This analysis indicated rpl32 is the only gene stable among the genes tested (SD value of 0.91). The comparative delta-Ct method was used to estimate the most stably expressed reference gene based on delta-Ct value variation. A lower value is considered more stable, and the results are similar to geNorm and BestKeeper with the gapdh, ubq, and ef1a being the most stable with a value of 1.34, 1.39, and 1.42. The final ranking suggests that the most stable reference gene across several larval tissues is gapdh followed by ubq, rpl32, ef1a, tubulin and sdf (Table 2).
For different adult tissues, the geNorm statistic algorithm indicated that gapdh and rpl32 are the most stable genes with an M score of 0.96, and tubulin and ef1a surpassed the cutoff value and are considered unstable, with M values of 1.6 and 1.79 respectively. The normfinder analysis also showed gapdh and rpl32 are also the most two stable genes, with a value of 0.54. For BestKeeper, only rpl32 did not surpass the cutoff value with an SD value of 1.04. The comparative delta-Ct method indicated that rpl32 and gapdh are the most stable genes, with a value of 1.46 and 1.48, respectively. The geometric mean ranking picked rpl32 and gapdh are the most stable genes across different adult tissues ( Table 2).
Summary of data from larval, adult (male and female) tissues and stages are included in Table 3. The geNorm statistic algorithm showed that gapdh and ubq are the most stable genes with an M score of 1.24, and ef1a. The gapdh gene is also identified as the most stable reference gene with a stability value of 0.74 and SD value of 1.44 using NormFinder and the comparative delta-CT algorithms, respectively. The BestKeeper method indicated that rpl32 and gapdh are the most stable genes, with a standard variation of 0.95 and 1.08. The final ranking calculated based on the combined algorithm values showed gapdh, rpl32, ubq, sdf, tubulin, and ef1a as the most to the least stable genes.
The gapdh gene identified as the most stable gene across different tissues of A. glabripennis larvae in the current study has also been identified as one of the most stable reference genes in other insects including Apis mellifera 33 , Spodoptera exigua and B. mori 22 . Since no stable reference gene has been identified in A. glabripennis, the reference genes identified in these studies would help A. glabripennis community with the gene expression studies.
Expression of iap. The expression of iap varied among the larval tissues tested. The highest mRNA levels were detected in the midgut followed by the foregut and hindgut. The lowest levels of IAP mRNA were detected in the integument followed by head and fat body (Fig. 3).

RNAi-response.
We prepared dsRNA targeting iap gene (dsIAP) and tested its efficiency in knocking down this gene in larvae. Five days after injection of dsIAP, the IAP mRNA levels were reduced by 81% when compared to their levels in control larvae injected with dsGFP (Fig. 4A). Injection of dsIAP also caused a significant mortality and only 10% of the treated larvae survived by 10 days after injection compared to 100% survival in control larvae injected with dsGFP (Fig. 4B). The knockdown of iap was also observed in specific tissues; fat body (77% knockdown; Fig. 5A) and alimentary canal (87% knockdown; Fig. 5B). RNAi response was also demonstrated in ALB adults. After dsIAP injection, a reduction of 91% in IAP mRNA levels was observed (Fig. 6A). The knockdown of iap gene in adults caused 100% mortality (Fig. 6B). Application of dsIAP caused mortality in adults and nymphs of Lygus lineolaris 34 .     Table 3. Primers used in identification of reference genes. R 2 : Correlation Coefficients; Eff: Amplification efficiency.
In conclusion, the identification of the core RNAi genes combined with RNAi bioassay results clearly demonstrates functional RNAi machinery in A. glabripennis. Based on our in silico findings, A. glabripennis possesses the core RNAi genes to successfully take up dsRNA and spread the signal within the insect. When dsRNA targeting an essential insect gene was injected into larvae and adults of A. glabripennis, nearly complete mortality of the treated insects occurred after target gene silencing in different insect tissues. These data further confirm the presence of functional RNAi machinery in A. glabripennis. Although the present study identified and validated the RNAi machinery in A. glabripennis, further studies to address the function of the AgladsRNase enzymes and their effects on RNAi response after ingestion of dsRNA need to be conducted. The characterization and understanding of the RNAi machinery could help in the development of RNAi-based methods to control this invasive   wood-boring insect pest. Also, the validation of reference genes will help with the normalization in future gene expression studies.

Material and Methods
Identification and phylogenetic analysis of RNAi genes. Core RNAi genes coding for proteins involved in miRNA, siRNA and piRNA pathways, including vesicle mediated transport, proton transport, intracellular transport, lipid metabolism and dsRNA uptake were identified in A. glabripennis. The amino acid sequences of the proteins were obtained from i5K workspace (http://i5k.nal.usd.gov/webapp/blast/) based on sequence homology searches by running BLASTp using the NCBI BLAST service (http://www.ncbi.nlm.nih. gov/) and UNIPROT BLAST service (http://www.uniprot.org/). These searches used T. castaneum sequences as the query. Hits (contigs) obtained from these search were used to verify their sequence similarity, identity, and detection of functional domains. Multiple sequence alignment and phylogenetic tree construction were performed using the sequences identified from T. castaneum (Coleoptera); Bombyx mori and Spodoptera litura (Lepidoptera); Apis mellifera (Hymenoptera); Droshophila melanogaster (Diptera); Acyrthosiphon pisum (Homoptera); Schistocerca americana and Locusta migratoria (Orthoptera); and the nematode Caenorhabditis elegans. The ClustalW program in MEGA 7.0 was used to align Dicer (Dcr1, Dcr2, Drosha), dsRBP (Pasha, R2D2, Loquacious), Argonaute (Ago1, Ago 2, Ago3, Aubergine), and Sid1-like protein (SilA, SilC) sequences. The neighbor-joining analysis was performed in MEGA 7.0 with bootstrapping to estimate the reliability of phylogenetic reconstruction (5000 replicates).
Insect rearing. A. glabripennis adults and larvae used in the studies were reared at the USDA-ARS Beneficial Insects Research Unit (BIIRU). The A. glabripennis colony was developed from beetles collected in 1999 from infested areas in New York, New Jersey, Chicago, and China and has since been maintained at 22.0 ± 2.5 °C, 45-75% RH, under a 14:10 h (L:D) photoperiod. Adult beetles were reared with fresh red maple twigs, and provided with larger logs for oviposition in 3.47-L plastic jars. To produce the appropriate instars of host larvae for tests, freshly cut red maple bolts were exposed to gravid A. glabripennis adults for one to two weeks in rearing jars (one female x one male per log per jar) with vented hard plastic lids under the rearing conditions mentioned above) to obtain appropriate levels of A. glabripennis egg deposition. Females of A. glabripennis chew oviposition pits and insert their eggs between the inner bark and the sap wood. The exposed logs were then incubated in the rearing room under the same environmental condition for three to four weeks until larvae hatched. Newly hatched larvae (approximately one to two weeks old after hatching) were then dissected out of infested logs and reared on artificial diet in 50 ml cup for four to six weeks before use in various experiments.
RNA extraction, primer design, and RT-qPCR. Total RNA was isolated from dissected tissue or whole larvae using the TRI Reagent RT (Molecular Research Center Inc.). The cDNA was synthesized using M-MLV Reverse Transcriptase (Invitrogen) from RNA and was used as a template in RT-qPCR. Candidate reference genes were selected based on the previous publications reporting stability in other insects 21,35,36 . Gene specific primers for each gene were designed using Primer3Plus (http://www.bioinformatics.nl/cgi-bin/primer3plus/ primer3plus.cgi). Sequences of primer are included in Table 4. PCR amplification efficiencies (E) and correlation coefficients (R 2 ) were checked to validate the RT-qPCR primers. Standard curves were constructed using 5-fold serially diluted cDNA, a mix of female and larva heads, for each primer pair. The expression analyses of the target genes were conducted using SYBR Green PCR Master Mix. Briefly, the PCR mixture contained 1 μL of synthesized cDNA, 0.2 μL of each primer (10 mM), 5 μL of the SYBR green PCR master mix and 3.6 μL of ddH 2 O. The reactions were carried out in triplicate per template in a final volume of 10 μL. RT-qPCR reactions were performed on the StepOnePlus Real-Time PCR system (Applied Biosystems) using the following cycling conditions: one cycle at 95 °C, followed by 40 cycles of denaturation at 95 °C, annealing and extension at 60 °C. At the end of each RT-qPCR reaction, a melting curve was generated to confirm a single peak and rule out the possibility of primer-dimer and non-specific product formation. For iap gene expression and gene silencing analyses, the (B) Survival rate (%) on 10 th day after injection of dsRNA. The asterisk indicates significant differences in mortality between treatment and control (Fisher's Exact test, Two-tailed, P-value < 0.001, N = 10).
2 −ΔΔCt method 37 was used to calculate the relative expression level of the target gene in the samples as compared to reference samples. For statistical analyses, one-way ANOVA, Tukey Test (P < 0.05) was used to test significance in differences of iap gene expression among different tissues.
Selection of reference genes. A web based tool, RefFinder (http://fulxie.0fees.us/?i=1) 26 , which integrates all four software algorithms, GeNorm 20 , NormFinder 23 , BestKeeper 24 and the comparative delta-Ct method 25 was used to evaluate reference gene stability from the experimental datasets. The mean Ct value of each sample and for each primer was used as the input data. dsRNA synthesis. Template DNA for dsRNA synthesis was amplified using the gene specific primers (Table 4). PCR conditions used were 94 °C for 4 min, followed by 35 cycles of 94 °C for 30 s, 60 °C for 30 s and 72 °C for 45 s, finishing with an extension step at 72 °C for 10 min. The PCR template was purified using a PCR purification kit (Qiagen). After PCR purification, dsRNA synthesis was performed using the MEGAscript RNAi Kit (Ambion Inc.) following manufacturer's instructions. Briefly, 200 ng of purified PCR product was used as template in a 20 μL in vitro transcription reaction. The reaction mix was incubated for 16 h at 37 °C, followed by 15 min of DNase treatment. The dsRNA was precipitated by adding 0.1 x volume of sodium acetate (3 M, pH 5.2) and 2.5 x volume of 100% ethanol and incubation at −20 °C for at least 2 h followed by centrifugation at 4 °C for 30 min. The dsRNA pellet was then rinsed with 750 μL of 75% ethanol and centrifuged again at 4 °C for 15 min. The ethanol was removed and the dsRNA was dried and diluted in ultrapure distilled water. The quality of the dsRNA was checked by electrophoresis and quantified using a spectrophotometer (NanoDrop Technologies). dsRNA targeting a fragment of the gene coding for green fluorescence protein was prepared using the protocol described above was used as a control. dsRNA bioassays. To assess gene silencing in A. glabripennis larvae, six to eight weeks-old larvae were injected with 2 μL containing 5 μg of dsRNA. For testing gene silencing in different tissues of A. glabripennis larvae (four to five weeks-old), and adults (two to three weeks after emergence) were injected with 2 to 2.5 μL containing 16 μg of dsRNA. Two treatments, dsIAP as a treatment and dsGFP as a control, were performed. After injection, adult beetles were maintained on freshly cut maple twigs and larvae were maintained on artificial diet in a plant growth chamber under normal rearing conditions described above. On the fifth day after dsRNA injection, the insects were collected and stored at −80 °C until RNA extraction. For statistical analysis, t-test [One-tailed (P ≤ 0.001)] was used for iap gene silencing compared to a single control (dsGFP injected). For the mortality bioassay, larvae (four to five weeks-old) and adults (one to two weeks-old) were injected with approximately 16 μg of dsRNA (in 2-2.5 μl) targeting iap or gfp gene. The mortality was scored on the 10 th day after the injection. For statistical analysis, Fisher's Exact test [Two-tailed (P < 0.001)] was performed.