Knockout of the HaREase Gene Improves the Stability of dsRNA and Increases the Sensitivity of Helicoverpa armigera to Bacillus thuringiensis Toxin

Double-stranded RNA (dsRNA)-induced genes are usually related to RNA interference (RNAi) mechanisms and are involved in immune-related pathways. In a previous study, we found a lepidopteran-specific nuclease gene REase that was up-regulated by dsRNA and that affected RNAi efficiency in Asian corn borer (Ostrinia furnacalis). In this study, to verify the function of REase, the homologous gene HaREase in cotton bollworm (Helicoverpa armigera) was knocked out using CRISPR/Cas9 system. We found that the midgut epithelium structure was apparently not affected in the ΔHaREase mutant [Knock out (KO)]. Transcript sequencing results showed that most of the known insect immune-related genes were up-regulated in KO. When second instar larvae were fed artificial diet with Cry1Ac, a protoxin from Bacillus thuringiensis (Bt), in sublethal doses (2.5 or 4 μg/g), the growth rate of KO was repressed significantly. The dsRNA stability was also enhanced in midgut extraction of KO; however, RNAi efficiency was not obviously improved compared with the wild type (WT). The KO and WT were injected with dsEGFP (Enhanced green fluorescent protein) and subjected to transcriptome sequencing. The results showed that the expression levels of 14 nuclease genes were enhanced in KO after the dsRNA treatment. These findings revealed that HaREase expression level was not only related with dsRNA stability, but also with Bt resistance in cotton bollworm. When HaREase was knocked out, other immune- or nuclease-related genes were enhanced significantly. These results remind us that insect immune system is complex and pest control for cotton bollworm is an arduous task.


INTRODUCTION
RNA interference (RNAi) technology is a potential strategy for crop protection against insect pests by double-stranded RNA (dsRNA) spraying or feeding, or via transgenic plants (Nandety et al., 2015;Joga et al., 2016;Zhang J. et al., 2017;Whitten, 2019). However, RNAi efficiency is relatively low in lepidopteran insects and this is a core factor that frustrates the application of this technology (Terenius et al., 2011;Cooper et al., 2019). Many factors affect RNAi efficiency in insects, including the absence of RNA-dependent RNA polymerase (RdRP)-mediated synthesis of secondary small interfering RNAs (siRNAs) (Sijen et al., 2001), the rate of dsRNA processing into siRNAs (Guan et al., 2018a), the dsRNA uptake and transport mechanism (Luo et al., 2012), and the degradation rate of dsRNA (Christiaens et al., 2014;Spit et al., 2017;Song et al., 2019).
In a previous study, we identified an RNAi efficiency-related nuclease gene REase that was induced specifically by dsRNA from Asian corn borer (Ostrinia furnacalis). The REase protein had nuclease activity and could degrade various types of nucleic acids in vivo and in vitro. Overexpression of REase decreased the RNAi efficiency in Drosophila, and knocking down the expression level of REase improved RNAi efficiency in Asian corn borer (Guan et al., 2018b). Therefore, the expression levels of dsRNAinduced genes can influence RNAi efficiency. Further studies of dsRNA-induced genes will help to better understand the RNAi and immune-related pathways, and guide the theory and practice of using RNAi in pest control.
In this study, to explore the function of REase in other insects, we cloned the homologous gene of HaREase in cotton bollworm by PCR. Then, we used the CRISPR/Cas9 system to knockout the HaREase in cotton bollworm. We performed transcriptome sequencing, dsRNA degradation tests, and Bacillus thuringiensis (Bt) sensitivity bioassays, to determine the effect of HaREase on RNAi-and immune-related pathways. The results revealed that HaREase was a lepidopteran-specific gene that was expressed mainly in the midgut and was induced by dsRNA. The HaREase expression levels affected dsRNA degradation through changed other nucleases. We also discovered that HaREase was involved in the lepidopteran immune stress processes and affected the resistance of cotton bollworm to Bt toxicity. These results will provide a novel strategy to enhance the sensitivity of insects to Bt toxin by inhibiting immune-related genes. These results will provide a novel strategy to enhance the stability of insects to Bt toxin by inhibiting immunerelated genes.

Insect Culture
Cotton bollworms were reared in the laboratory at 25 • C and 75% relative humidity on a 14-h/10-h light/dark cycle. The larvae were fed a modified artificial diet (120 g of maize granules, 32 g of maize flour, 120 g of soybean flour, 4 g of vitamin C, 12 g of agar, 72 g of yeast powder, 4 g of sorbic acid, 60 g of glucose, 1.6 ml of formaldehyde, and 1000 ml of water). Moths were fed a 10% (v/v) honey solution.

Single Guide RNA (sgRNA) Design and Synthesis
Ten individual cotton bollworms were used to identify of the conservative DNA fragment of the target HaREase gene by PCR with primers pairs (5 -GCCAAAGAAGAATCAAAGACTTATC AT-3 and 5 -CCTCCTTGGCTTCAGTGACGCAATACC-3 ). One small conservative sequence 5 -GCTGATGAAAGAT AGTCTCGCGG-3 was selected as the target site and the corresponding sgRNA was synthesized using GeneArt TM Precision gRNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, United States). Basically, the target sequence of 5 -N20-3 before PAM sequence "CGG" was used to design the sgRNA. A pair of long primers (5 -TAATACGACTCACTATAGCTGATGAAAGATAGT-3 and 5 -TTCTAGCTCTAAAACCGAGACTATCTTTCATCA-3 ) was designed and the sgRNA was synthesized following the manufacturer's instructions. The Cas9 protein (GeneArt TM Platinum TM Cas9 Nuclease) was obtained from Thermo Fisher Scientific.

Microinjection and Screening of Mutants
Within half an hour of oviposition, the freshly laid eggs were collected and placed on a microscope slide and fixed with double-sided adhesive tape. Cas9 protein and sgRNA with final concentrations of 200 and 500 ng/µL, respectively, were microinjection into these pretreated eggs. Approximately 600 eggs were injected and incubated at 25 • C and 75% relative humidity for 3-4 days until hatching. When these injected F0 individuals grew into adults, males and females were single paired with wildtype individuals in plastic cups for mating and laying eggs. The newly hatched F1 larvae were maintained into adults and single paired by inbreeding within same family. The laidegg-adults of F1 were collected and their whole bodies were used for mutation screening by sequencing the genome DNA of each candidate individuals. If both of the F1 adults contain same mutant genotype, homozygous of mutants will be obtained in the F2 generation. The HaREase KO lines for functional analysis were obtained by inbreeding using F2 homozygous mutants. The genomic DNA was isolated using a TIANamp Genomic DNA Kit (Tiangen, Beijing, China). PCR products were purified using an EasyPure PCR Purification Kit (TransGen Biotech, Beijing, China), cloned into a pEASY-Blunt cloning vector (TransGen Biotech, Beijing, China), and sequenced for mutation screening.
dsRNA Preparation dsRNAs were synthesized using a MEGAscript RNAi kit (Ambion, Huntingdon, United Kingdom) in accordance with the manufacturer's instructions. T7 promoter sequences were tailed to the 5 -ends of the DNA templates by PCR amplifications. All the primer sequences are listed in Supplementary  Table S1. Template DNA and single-stranded RNA was removed from the transcription reaction by DNase and RNase treatments, respectively.
In vitro Degradation of HaREase on dsRNA Midgut extraction samples (Guan et al., 2018b) were collected from the fourth instar wild type (WT) and HaREase mutant (KO), put them into 100 µl PBS, puncture and release the midgut fluid, take the supernatant after centrifugation. Dilute the midgut fluid with a final protein concentration of 3 ng/µl. After that, dsRNA was incubated with the diluted midgut fluid samples with a final concentration of 50 ng/µl at 37 • C. The dsRNA samples were collected at different time points and the dsRNA integrity was analyzed by 1% agarose gel electrophoresis. The experiment was repeated four times.

Sample Collection and RNA Isolation
Cotton bollworm samples were collected, immediately frozen in liquid nitrogen, and stored at −80 • C until RNA extraction. Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) in accordance with the manufacturer's instructions. Samples were treated with RNase-free DNaseI (New England Biolabs, Ipswich, MA, United States) for 30 min at 37 • C to remove residual DNA.

Real-Time Quantitative PCR (RT-qPCR)
Total RNA was extracted using TRIzol R reagent (Invitrogen, Carlsbad, CA, United States) in accordance with the manufacturer's instructions. First-strand cDNA was synthesized from 1 µg of RNA primed by oligo(dT) 18 using M-MLV reverse transcriptase (Takara, Kyoto, Japan). The RT-qPCR assay for multiple genes was performed using SYBR R Premix Ex TaqTM II (Takara, Kyoto, Japan). The primer sequences are listed in Supplementary Table S1. Melting curve analyses were performed for all of the primers. 18S rRNA expression levels were used to normalize the cycle threshold (Ct) (Chandra et al., 2014). The RT-qPCR was carried out using a Mastercycler ep realplex (Eppendorf, Hamburg, Germany). All the RT-qPCR assays were repeated three times. To assess the extent of RNAi, RNA was extracted from pools of three dsRNA-treated surviving larvae using TRIzol reagent (Invitrogen, Carlsbad, CA, United States), and each treatment was repeated three times. RT-qPCR reactions and data were analyzed in accordance with the methods of Livak and Schmittgen (2001). The data were analyzed using one-way analysis of variance to detect treatment effects compared with the untreated control.

Transcriptome Sequencing and Analysis
The third instar larvae of WT and KO cotton bollworm were injected with dsEGFP (10 µg/larvae). Six hours later, samples were collected and frozen in liquid nitrogen for RNA sequencing (RNA-seq). Each sample contained 10 larvae. The non-treated samples were used as the control (WT). All the samples were sequenced using an Illumina HiSeq 2000 analyzer at BGI, Shenzhen, China. The datasets generated during the current study have been deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA548442.
Sequences were searched against the GenBank non-redundant database (Nr) using the BLASTx algorithm. Gene ontology (GO) annotation was assigned using default annotation parameters in Blast2go 1 . The Cluster of Orthologous Groups (COG) and Kyoto encyclopedia of Genes and Genomes (KEGG) pathway annotations were assigned using the BLASTall online program against the respective databases 2,3 . The RNA-seq expression levels were calculated using the reads per kilobase per million mapped read (RPKM) method.

Light Microscopy
Tissue sectioning was performed as described by Fischer et al. (2008). Slides (5 µm) were stained using hematoxylin and eosin solution to show morphological changes. All images were taken using a fluorescent microscope (BX51, Olympus, Tokoyo, Japan) with differential interference contrast and the appropriate filter, unless otherwise stated.

Bt Toxin and Bioassays
The Cry1AC protoxin was obtained from JZ Wei at Henan Agricultural University. Diet overlay bioassays (Hernandez-Rodriguez et al., 2008) were used to evaluate larval response to the Bt toxins. Cry1AC was diluted in 50 mM Na 2 CO 3 (pH 10.0) and then added to the diet to give final concentrations of 2.5 µg or 4 µg Cry1AC/g diet. We used 12-well plates for the bioassay. Second instar larvae that hatched within 6 h were transferred onto the diet surface in each well. Each treatment was performed in triplicate (12 larvae in each plate was regarded as one repeat). After 3 days, the larvae were weighed and the statistical analysis was performed.

Spatiotemporal Expression Patterns of HaREase in Cotton Bollworm
The full-length 1737-nucleotide cDNA sequence of HaREase was cloned by PCR (XM_021337058.1). Then, the spatiotemporal expression levels of HaREase at eight development stages from egg to adult stages, and in 10 different tissues (head, epidermis, foregut, midgut, hindgut, hemolymph, fat body, malpighian tube, salivary gland, ganglion) were investigated by RT-qPCR (Figures 1A,B). The results indicated that HaREase was expressed mainly in the guts of larvae from the third to fifth instar larval stages. The expression pattern of HaREase in cotton bollworm was similar to that of REase in Asian corn borer. HaREase also was up-regulated by dsEGFP just as REase was in Asian corn borer ( Figure 1C; Guan et al., 2018b). These results imply that REase genes may have similar functions in different insect species.

Generation of HaREase Knockout Mutants of Cotton Bollworm
The CRISPR/Cas9 system was used to knock out HaREase for further function analysis. sgRNA target site that covered the HaREase coding region was selected following the 5 -N20-NGG-3 rule (Figure 2A). For the screening of mutants, six mutant lines with different genotypes were obtained in which HaREase-1, HaREase-2, HaREase-4, and HaREase-6 were effective mutants, but HaREase-3 and HaREase-5 were not because they did not break the ORF coding for the proteins. At last, the HaREase-6 mutant line which lack two bases of nucleic acid sequences was selected for all the subsequent experiments ( Figure 2B). Gene silencing of HaREase-6 KO mutants was confirmed by RT-qPCR, which showed that the mRNA level of HaREase was significantly lower than in WT ( Figure 2C). These results indicated that HaREase was effectively knocked out in this strian.

Transcriptome and Intestinal Microorganism Analysis Between KO and WT Cotton Bollworm
We noticed that KO had no obvious visible differences in insect growth and development compared with WT. Because HaREase was expressed mainly in the gut of larval stage cotton bollworm (Figures 1A,B), we dissected out the midgut of fifth instar larvae of KO and WT individuals to look for anatomical differences. No obvious differences were observed ( Figure 2E). Furthermore, analysis of the 16s rDNA from the intestinal microorganism of third instar larvae by RT-qPCR found no significant effect on intestinal flora of KO compared with WT ( Figure 2D). To investigate the effect of HaREase knock out on other genes, we performed transcriptome sequencing analysis of the third instar larvae midguts of WT and KO. Compared with WT, 3574 transcripts were up-regulated in KO (LogFC > 1, FDR < 0.005) and 325 transcripts were down-regulated (LogFC < −1, FDR < 0.005) (Figures 3A,B). The KEGG pathway analysis of these differentially expressed genes identified spliceosome, protein processing in endoplasmic reticulum, ubiquitin mediated proteolysis, and mRNA surveillance as the top metabolic pathways associated with the knockout of HaREase (Figure 3C). Four immune-related pathways were further calculated, unexpectedly, most of the immune-related genes were up-regulated in KO strain ( Table 1).

HaREase Knockout Enhanced the Sensitivity of Cotton Bollworm to Bt Toxin
To understand the effects of the up-regulated immune-related genes on cotton bollworm, the larvae were treated with sublethal doses of Bt toxin. Second instar larvae were fed an artificial diet containing the Cry1Ac protoxin in sublethal doses of 2.5 or 4 µg/g. Three days later, the larvae were weighed, and the weight gain of the KO larvae was significantly lower than that of the WT larvae (Figures 4A,B). This indicated that HaREase may be involved in insect immunity or detoxification of Cry1Ac.
After feeding with Cry1Ac, the expression levels of several immune related genes were detected, cecropin-1, galiomicin, gloverin, lysozyme, ceropin-3, cobatoxin, and moricin were upregulated significantly in both WT and KO (Figure 4C), whereas GNRP and cecropin-2 were up-regulated only in WT. This may be related to difference in Cry1Ac resistance between KO and WT strain ( Figure 4C).
By comparing the expression levels of Bt receptor and resistance-related genes between WT and KO, we detected one Bt receptor gene and three genes encoding resistance-related enzymes that were down-regulated significantly in KO ( Table 2). Three other Bt receptor genes were also up-regulated significantly in KO, especially the major Bt receptor gene that encodes aminopeptidase N ( Table 2). These results may explain why KO was more sensitive to Bt toxin than WT, and indicate HaREase as a new candidate gene for tackling the Bt resistance in cotton bollworm.

Effect of HaREase on dsRNA Stability and the RNAi Pathway
To further confirm the effect of HaREase on dsRNA stability, the midgut fluid was extracted from the fouth instar larvae of WT and KO. dsRNA was incubated in the midgut extract and agarose gel electrophoresis was used to detect the degradation rate of the dsRNA. The result showed that dsRNA degraded more quickly in the midgut extract from WT compared with the midgut extract from KO ( Figure 5A). NUDFV2 was selected to test RNAi efficiency between WT and KO. The same amount of dsNUDFV2 was injected into third instar larvae of WT and KO. After 6 h, samples were collected and the expression level of NUDFV2 was checked. However, the RNAi efficiency of dsNUDFV2 was not enhanced obviously in KO ( Figure 5B). Next, we checked the expression levels of the core RNAi pathway genes argonaute-2 and dicer-2 by RT-qPCR and found that both these genes were up-regulated in KO after dsRNA treatment, but no significant differences in their expression levels were detected between WT and KO after the dsRNA treatments (Figures 5C,D). These results indicated that HaREase knockout had no obvious effect on RNAi efficiency in cotton bollworm, and are completely different from what was found in Asian corn borer (Guan et al., 2018b). This suggests that the RNAi mechanism in cotton bollworm as a  polyphagous insect is more complicated and variable than in Asian corn borer.
To further clarify the relationship between HaREase and RNAi efficiency, the third instar cotton bollworm larvae of WT and KO were injected with dsEGFP (10 µg/larvae). After 6 h, the samples were collected for transcriptome sequencing. The results indicated that many genes had changed after dsEGFP treatment, among them, 1176 transcripts were co-changed genes ( Figure 5E). The KEGG annotations indicated that these genes were involved mainly in spliceosome, protein processing in endoplasmic reticulum, and ubiquitin mediated proteolysis pathways ( Figure 5F). The heat map showed that the expression levels of some of the transcripts were significantly different between the WT and KO samples ( Figure 5G).
We thoroughly analyzed the nuclease genes in WT and KO and found that 12 of 15 nuclease genes were up-regulated in KO before and after dsRNA treatment compared with their expression in WT ( Table 3). The increased expression of these nucleases in KO may explain why the RNAi efficiency of cotton bollworm was not enhanced by knocking out HaREase.

DISCUSSION
The KO mutants that we obtained using CRISPR/Cas9 technology were used to study the function of HaREase in cotton bollworm. The results indicate that HaREase is a lepidopteranspecific gene that can be induced specifically by dsRNA and is FIGURE 4 | Sensitivity of KO and WT cotton bollworm to Bt toxins. (A,B) Growth and development of KO and WT cotton bollworm after feeding with an artificial diet containing Cry1Ac. (C) Expression levels of immune-related genes in KO and WT after feeding with Cry1Ac. Data are mean SD, n = 3. * P < 0.05; * * P < 0.01.  expressed mainly in the gut of cotton bollworm at the larval stage. The HaREase mutation had no obvious effect on the growth and development of cotton bollworm, but KO was more sensitive to Cry1Ac than WT. The results suggest that HaREase may be associated with insect immunity and will be a promising new target for tackling the Bt resistance in cotton bollworm. A primary concern is the impact of HaREase mutations on the RNAi efficiency of cotton bollworm. Many results have indicated that the nucleases used for dsRNA degradation can affect RNAi efficiency. Knock-down of dsRNase in the Leptinotarsa decemlineata midgut enhanced RNAi efficiency (Spit et al., 2017). Our previous research indicated that the expression level of REase was related to RNAi efficiency in Asian corn borer (Guan et al., 2018b). Therefore, we tested the dsRNA degradation rate in midgut extractions from WT and KO and found that although dsRNA stability was enhanced in KO, the RNAi efficiency was not improved obviously (Figure 5).
RNA interference efficiency depends on the dsRNA degradation rate by insect nucleases (Wang et al., 2016;Song et al., 2019). Nuclease activity has been detected in Culex pipiens quinquefasciatus (Calvo and Ribeiro, 2006), Glossina morsitans (Caljon et al., 2012), Nezara viridula (Lomate and Bonning, 2016), Bombyx mori (Arimatsu et al., 2007), and Schistocerca gregaria (Wynant et al., 2014). Many researches have revealed that RNAi efficiency in different insect orders species was highly correlated with dsRNA stability (Wang et al., 2016;Song et al., 2019). The dsRNA was more stable in coleopteran insects than in lepidopteran insects and RNAi efficiency was higher in coleopteran than in lepidopteran (Ivashuta et al., 2015). Therefore, knocking down nuclease activity using RNAi technology can enhance RNAi efficiency (Spit et al., 2017;Guan et al., 2018b;Song et al., 2019). For instance, dsRNase knockdown enhanced RNAi efficiency in L. decemlineata (Spit et al., 2017). However, this does not work in all insect species, and dsRNase knock down did not enhance RNAi efficiency in S. gregaria (Spit et al., 2017). The diversity of nucleases in different insects may be the main reason for the observed differences in RNAi efficiency. In this study, the transcriptome sequencing results confirmed that when one nuclease gene HaREase was knocked out, the expression levels of a large number of other nuclease genes were enhanced to compensate. Therefore, at least for cotton bollworm, it is difficult to increase RNAi efficiency by knocking out only one nuclease gene. When exogenous dsRNA is introduced into the insect body, it can be recognized as a pathogen-associated molecular pattern (PAMP) by pattern-recognition receptors (PRRs) to induce an immune response (Sadler and Williams, 2008). Indeed, the RNAi pathway is a natural antiviral immunity system (Ding, 2010;Gammon and Mello, 2015;Meki et al., 2018). The RNAi pathway related genes argonaute-2, dicer-2 and the immunerelated genes hemolin, attacin, gloverin, lysozyme, and NF-kappa B are induced by dsRNA (Hirai et al., 2004;Guan et al., 2018c). Because HaREase is expressed mainly in insect gut and can be induced by dsRNA, we speculated that HaREase may be related to the insect immune response. Knockdown of an immune-related gene in Spodoptora littoralis increased its susceptibility to Bt toxin (Caccia et al., 2016). When we combined Bt toxin with dsRNA targeting chymotrypsinlike genes, larval mortality increased and inset death was accelerated (Guan et al., 2016). Thus, a new strategy for enhancing the impact of Bt toxins on insects may be the knocking down of immune-related genes. When sublethal doses of Cry1Ac were added into the artificial diet, the KO mutants were more sensitive to Bt and their weight gains decreased significantly (Figures 4A,B), indicating the expression level of HaREase was related to insect sensitivity to Bt toxin. Therefore, HaREase may be a new target for improving Bt sensitivity in cotton bollworm.
Bt toxins from B. thuringiensis have been used as an insecticide in cotton bollworm, which has mitigated the impact of the destructive cotton bollworm pest worldwide (Sanahuja et al., 2011;Romeis et al., 2019). When Bt protoxins are eaten by insects, they can be digested by midgut proteases into activated toxins that bind to receptors in the insect midgut, forming lytic pores in the membrane and leading to cell breakdown and insect death (Gill et al., 1992;Bravo et al., 2011). However, the evolution of insect resistance to Bt is gradually reducing the benefits of this approach (Tabashnik and Carriere, 2017). Many factors are related to Bt sensitivity, including the Bt receptors and the proteases involved in protoxin activation (Rajagopal et al., 2009;Cao et al., 2013;Liu et al., 2014;Jin et al., 2018). For example, the reduced binding capacity of Bt toxins to midgut receptors such as like ALP, APN, cadherin, and ABCC2 (Xu et al., 2005;Jurat-Fuentes and Adang, 2007;Gahan et al., 2010;Baxter et al., 2011;Atsumi et al., 2012;Xiao et al., 2014;Chen et al., 2015;Zhang T. et al., 2017), and the reduced expression levels of trypsin genes, and ALP2 are all related to Bt resistance in cotton bollworm (Wei et al., 2018). In this study, we identified three Bt receptor genes that were enhanced significantly in KO ( Table 2). Besides this, after feeding with a diet containing Cry1Ac, the expression levels of the immune-related genes GNRP, PGRP-C, cecropin-1, galiomicin, cecropin-2, cobatoxin, moricin, and gallerimycin were lower in KO than in WT ( Figure 4C). Together, these results suggest that HaREase is involved in insect immune stress response against Bt toxicity.
The insect immune response against various external invasions is a complex process, and some immune genes are lepidopteran-specific. Hemolin is a lepidopteran-specific gene that encodes a pattern recognition protein that can be induced by bacteria and viruses (Aathmanathan et al., 2018). HaREase is also lepidopteran-specific gene (Guan et al., 2018b), study these lepidopteran specific genes will help in understanding the specific immune pathway of these insects and contribute to the management of cotton bollworm.

DATA AVAILABILITY STATEMENT
The datasets generated during the study have been deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA548442.

AUTHOR CONTRIBUTIONS
RG and QC performed the experiments. RG and BY designed and wrote the manuscript. HL and SH analyzed the data.
XM and GW revised the manuscript. All authors read and approved the manuscript.

FUNDING
This work was supported by the National Natural Science Foundation of China (31702057, 31672354, and 31772520), and the National Key R&D Program of China (2017YFD0200900). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.