The Myotis chiloensis Guano Virome: Viral Nucleic Acid Enrichments for High-Resolution Virome Elucidation and Full Alphacoronavirus Genome Assembly

Bats are widespread mammals of the order Chiroptera. They are key for ecosystem functioning, participating in crucial processes. Their unique ability amongst mammals to fly long distances, their frequently large population sizes, and their longevity favor infectious agent persistence and spread. This includes a large variety of viruses, encompassing many important zoonotic ones that cause severe diseases in humans and domestic animals. Despite this, the understanding of the viral ecological diversity residing in bat populations remains unclear, which complicates the determination of the origins of zoonotic viruses. To gain knowledge on the viral community of a widely distributed insectivorous bat species, we characterized the guano virome of a native Chilean bat species (Myotis chiloensis (Waterhouse, 1840)). By applying a novel enrichment strategy, we were able to secure a consequent percentage of viral reads, providing unprecedented resolution for a bat virome. This in turn enabled us to identify and assemble a new bat alphacoronavirus from Chilean bats closely related to PEDV, an important viral pathogen with high mortality rates in suckling piglets. This study highlights the importance of applying and improving high-resolution virome studies in this vital order to ultimately enhance epidemiological surveillance for potentially zoonotic pathogens.

Introduction
Bats are mammals in the order Chiroptera that can be found in all continents except the polar regions and some oceanic islands [1]. They are the second most widely distributed mammalian group after rodents and are represented by at least 1400 extant species worldwide [1][2][3]. They are a key group for ecosystem functioning, participating in crucial processes, including pollination and seed dispersal; control through predation of vertebrate and invertebrate populations (comprising pest species), indirectly controlling the composition and structure of plant communities through its effect on herbivory; and recycling of energy, being able to mobilize large volumes of nutrients [4,5]. Additionally, bats are hosts of a wide diversity of parasites and pathogens, and their diversity can have profound implications for the transmission of diseases in ecosystems.
Infectious agents exist in commensal, mutual, or parasitic relationships with their hosts. The relationship is restricted to a limited number of host species due to genetic adaptations developed during co-evolution. Spillage of pathogens between different animal species can occur sporadically. This can result from accidental or deliberate intrusion by one animal species into the ecological niche of another [6].
Bats are unique among mammals in their ability to fly and inhabit diverse ecological niches. Their ability to fly long distances paired with their longevity favor infectious agent persistence and spread [7]. In addition to this, changing environments highlight the shared health risks by humans and bats. These characteristics together with their frequently large populations emphasize their potential as hosts for pathogens [3]. Their role in disease epidemiology is further supported by the diversity of microorganisms detected in their microbiomes such as bacteria, fungi, parasites, and viruses. So far, viruses belonging to at least 14 families have been identified in bats [8]. Previous and ongoing research is predominantly focused on viral agents. This is mostly due to the potential for cross-species contamination, as illustrated by the Ebola outbreak in West Africa from the mid 2000 s and the COVID-19 pandemic [9]. Identification in bats of additional relevant viral pathogens with cross-species transmission potential such as rabies virus [3,10], henipaviruses (Nipah and Hendra viruses), and coronaviruses (CoVs) has further highlighted the importance of monitoring and characterizing the viromes of these key players in a given ecosystem [11][12][13].
CoVs are classified in four different genera: Alphacoronavirus, Betacoronavirus, Gammacoronavirus, and Deltacoronavirus [14]. These viruses can infect humans, domestic and wild mammals and birds, causing respiratory and enteric diseases. Surveillance programs around the world have made an effort to better understand bats' role as reservoirs for several viruses that pose significant threats to human health, the poultry industry, and hog farming, amongst others. Several different bat species share an ecosystem with humans, for example, Myotis chiloensis [15], Tadarida brasiliensis [16], and farm animal species, which can favor zoonotic infections [17]. Swine Acute Diarrhea Syndrome Coronavirus (SADS-CoV) and porcine epidemic diarrhea virus (PEDV) are examples of coronaviruses that jumped into piglets, causing acute gastroenteritis in neonatal piglets and ultimately, important economic losses [18]. In silico analyses suggest that SADS-CoV and PEDV may have originated from bat CoVs, reaffirming interspecies transmission and further suggesting the importance of viral surveillance programs in bats [11].
Myotis chiloensis is a small Vespertilionid bat species endemic to the southern cone of South America. M. chiloensis feeds on insects-primarily nematocerans-that are captured in flight [19,20]. They are currently considered a "Least Concern" species by the International Union for Conservation of Nature [19]. M. chiloensis is usually found cohabitating with other bat species, in particular with the Brazilian free-tailed bat (Tadarida brasiliensis), the big-eared brown bat (Histiotus macrotus), the small big-eared brown bat (H. montanus), and the southern big-eared brown bat (H. magellanicus) [20][21][22]. Although some studies have described the bacterial and fungal [2,23] communities that they harbor, little is known about the virome of these widely distributed insectivorous bats.
In order to gain knowledge on the viral community in a widely distributed insectivorous bat species, we characterized the guano virome of M. chiloensis. We monitored a bat colony located in the O'Higgins region of Chile from 2018-2021 ( Figure 1). By applying a combined enrichment strategy, we were able to secure a consequent percentage of viral reads. This entailed pairing flocculation for the enrichment of viral particle concentration, molecular enrichment through automated depletion, and bioinformatic filtering of contaminant sequences. Emerging infectious diseases of zoonotic origin embody a major public health and animal industry challenge. Here we describe a new bat alphacoronavirus from a Chilean bat closely related to PEDV an important viral pathogen with high mortality rates in suckling piglets, whose zoonotic potential should be further analyzed. We also highlight the importance of maintaining viral surveillance in this field.

Material and Methods
Sample collection: Fecal samples were obtained from bat colonies established in human buildings between the Metropolitan Region (MR) and O'Higgins Region (VI) of Chile at different dates between 2019 and 2021 ( Figure 1). They were stored at −80 °C. It should be noted that both samples from location C collected in 2019 were pooled for sequencing in order to get an overview of the colony's virome.
PCR screening: Fecal samples were homogenized in viral transport medium (VTM) and centrifuged at 6000× g for 15 min at 4 °C, and supernatant was stored at −80 °C. Viral RNA was extracted from 250 µL of the supernatant obtained with TRIzol ® Reagent (Thermo Fisher Scientific Waltham, MA, United States) based on the method of Chomczynski and Sacchi [24]. The RNA obtained was stored at −20 °C. Coronavirus screening was performed by a pancoronavirus one-step RT-PCR based on primers for a 251 bp conserved region in the polymerase gene (Cor-FW 5′-ACWCARHT-VAAYYTNAARTAYGC−3′; Cor-RV 5′-TCRCAYTTDGGRTARTCCCA−3′) [25]. RT-PCR was made with RT-PCR KAPA SYBR ® FAST one-step kit (KAPA Biosystems, Wilmington, MA, United States). [25]. The PCR product was run on a 1.5% agarose gel in TAE 1x at a constant 80 V for 30 min. Paramyxovirus screening was made by a panparamyxovirus one-step RT-PCR based on primers for a 121 bp conserved region in domain III of the RNA-dependent RNA polymerase gene (PMX1-FW 5′-GAR-GGIYIITGYCARAARNTNTGGAC−3′; PMX2-RV 5′-TIAYIGCWATIRIYTGRT-TRTCNCC−3′) [26]. RT-PCR was performed with Brilliant III Ultra-Fast qRT-PCR Master Mix (Agilent Technologies, Santa Clara, CA, United States). PCR product was run on a 2.5% agarose gel in TAE1x at a constant 80 V for 30 min.
Viral enrichment: We weighed 1 g of bat guano and resuspended it until reaching 5 mL of PBS (the ratio was 1:5, buffer was previously autoclaved and filtered using 0.2 µM esterase cellulose filters). Samples were vortexed vigorously, then centrifuged at 2000× g for 2 min at 4 °C, taking the supernatant. Next, the supernatant was centrifuged at 400× g for 10 min at 4 °C; taking once again the supernatant, this step was repeated once [27]. Next, the sample was centrifuged at 16,000× g for 45 min and the supernatant was kept. Samples were then filtered using a 0.2 µM sterivex made out of PES using a 50-milliliter syringe. Samples were then treated with DNase (Goldenbio, Montclair, CA, United States,

Material and Methods
Sample collection: Fecal samples were obtained from bat colonies established in human buildings between the Metropolitan Region (MR) and O'Higgins Region (VI) of Chile at different dates between 2019 and 2021 ( Figure 1). They were stored at −80 • C. It should be noted that both samples from location C collected in 2019 were pooled for sequencing in order to get an overview of the colony's virome.
PCR screening: Fecal samples were homogenized in viral transport medium (VTM) and centrifuged at 6000× g for 15 min at 4 • C, and supernatant was stored at −80 • C. Viral RNA was extracted from 250 µL of the supernatant obtained with TRIzol ® Reagent (Thermo Fisher Scientific Waltham, MA, United States) based on the method of Chomczynski and Sacchi [24]. The RNA obtained was stored at −20 • C. Coronavirus screening was performed by a pancoronavirus one-step RT-PCR based on primers for a 251 bp conserved region in the polymerase gene (Cor-FW 5 -ACWCARHTVAAYYTNAARTAYGC−3 ; Cor-RV 5 -TCRCAYTTDGGRTARTCCCA−3 ) [25]. RT-PCR was made with RT-PCR KAPA SYBR ® FAST one-step kit (KAPA Biosystems, Wilmington, MA, United States). [25]. The PCR product was run on a 1.5% agarose gel in TAE 1x at a constant 80 V for 30 min. Paramyxovirus screening was made by a panparamyxovirus one-step RT-PCR based on primers for a 121 bp conserved region in domain III of the RNA-dependent RNA polymerase gene (PMX1-FW 5 -GARGGIYIITGYCARAARNTNTGGAC−3 ; PMX2-RV 5 -TIAYIGCWATIRIYTGRTTRTCNCC−3 ) [26]. RT-PCR was performed with Brilliant III Ultra-Fast qRT-PCR Master Mix (Agilent Technologies, Santa Clara, CA, United States). PCR product was run on a 2.5% agarose gel in TAE1x at a constant 80 V for 30 min.
Viral enrichment: We weighed 1 g of bat guano and resuspended it until reaching 5 mL of PBS (the ratio was 1:5, buffer was previously autoclaved and filtered using 0.2 µM esterase cellulose filters). Samples were vortexed vigorously, then centrifuged at 2000× g for 2 min at 4 • C, taking the supernatant. Next, the supernatant was centrifuged at 400× g for 10 min at 4 • C; taking once again the supernatant, this step was repeated once [27]. Next, the sample was centrifuged at 16,000× g for 45 min and the supernatant was kept. Samples were then filtered using a 0.2 µM sterivex made out of PES using a 50-milliliter syringe. Samples were then treated with DNase (Goldenbio, Montclair, CA, United States, 100 U/mL of the original sample (from step 1), 100 U/mL RNAse (Goldenbio) A, and 100 U/mL Turbonuclease (Sigma, San Luis, MO, United States). Next, samples were incubated at 37 • C for 1.5 h or at 4 • C O/N. RNA shield was then added for enzyme inactivation. Samples were concentrated using PEG-NaCl and spinning at 4000× g for 30 min and then through ultracentrifugation at 100,000× g for 1 h. Total RNA was then extracted using the QIamp viral RNA extraction kit. Total RNA QC check and quantification using the nanochips for the bioanalyzer and finally, concentration was measured using Invitrogen's Qubit.
Depletion and library preparation and sequencing: After total RNA enrichment for viral sequences through centrifugation, total RNA was depleted for rRNA using Illumina's Ribozero Plus kit with a custom protocol in the Magelia (Inorevia, Paris, France). After bioinformatic verification of probe compatibility for bacterial/mammalian selection, 2 µLs of total RNA were used in the instrument for fully automated molecular enrichment of viral sequences. After depletion, total RNA was run in the Bioanalyzer 2100 using a nano RNA chip to determine rRNA contamination. Next, the total stranded RNA illumina protocol for library preparation was performed using 100 ng of depleted RNA. Library preparation yielded 98.2 nM of high-quality sequenceable material that was then run in a NextSeq 500 Mid Output Kit v2 (150 cycles) cartridge. This yielded a total of 90,530,942 raw reads. The data for this RNA-seq study has been submitted in the GenBank Sequence Read Archive (SRA) under project number PRJNA781397, SUB10683119.
Bioinformatic analyses: Reads were trimmed and read quality was assessed using the QC for sequencing reads tool from the Genomics Workbench (CLC, Qiagen, Paris, France). Mapping and de novo assembly were also performed using this software. Reads were then mapped to three different Myotis spp. genome assemblies in order to remove sequences of bat origin. Then, de novo assembly was performed on high quality reads under the following parameters: using updated contigs with an automatic bubble size, minimum contig length = 200, applying automatic word size and scaffolding, the paired distances were automatically detected, the mismatch cost was 2, insertion cost 3, deletion cost 3, length fraction 0.5, and similarity fraction 0.5. Finally, reads were mapped back to the contigs to consolidate them.
The assembled contigs were then automatically blasted under the following parameters using the blast nr database at NCBI: Match/Mismatch and Gap Costs = (Match 2 Mismatch 3 Existence 5 Extension 2), Expectation value = 0.05, Word size = 11, masking low complexity regions, a maximum number of hits at 100 and opening the search to all organisms. Assembly of the viral sequences was additionally manually curated after determining ORFs.
In parallel, sequences were treated using the Kaiju web server for taxonomic classification under the following parameters: we used the NCBInr protein database, low complexity protein sequences were filtered, Geedy and MEM run modes were both applied, minimum match length was 11, and minimum match score was 75 with 5 mismatches allowed.
Phylogenetic analysis: To perform phylogenetic analyses, the full-length sequence of the newly assembled M. chiloensis alphacoronavirus 1 was aligned with assorted related viral sequences isolated in different species including bat, murine, feline, porcine, avian, civet, and human origin. A total of 43 sequences were used. The accession numbers can be found in Table S1. The algorithm used for the alignment is the QIAGEN Aarhus applying the "very accurate" setting. The alignment was performed under the following parameters: Gap open cost = 10, Gap extension cost = 1,0, End gap cost = As any other. This alignment was then used in the IQ-TREE webserver to generate a phylogenetic tree. Modelfinder was paired with tree reconstruction and ultrafast bootstrap (1000 and 10,000 replicates). The Best-fit model according to BIC was GTR + F + I + G4. The model of rate heterogeneity was Invar + Gamma with 4 categories. Maximum likelihood and consensus trees were generated. Full results for this analysis can be found in Supplementary file 1, including bootstrap values.

Collected Samples
Samples were collected from four different locations from two different regions in central Chile. Locations A, B, and D were negative for viral screening; however, location C was positive for paramyxovirus for samples collected in 2018 and 2021. In addition to this, two samples from 2019 were positive for CoV, and these samples were pooled, purified, and ultraconcentrated for molecular enrichment through depletion and NGS library prep/sequencing (Figure 1).

Sample Enrichment and Virome Elucidation/Alphacoronavirus Genome Assembly and Phylogeny
Depletion of previously enriched total RNA extracted from guano samples allowed us to reach a 0% rRNA contamination level after Bioanalyzer verification using a Nanochip. After successful library prep and sequencing, a total of 90,530,942 high-quality reads were obtained. These reads were filtered through three different Myotis spp. assemblies prior to de novo assembly of these reads. This enabled us to generate a total of 5460 contigs. These contigs were blasted in NCBI's nr database, significant identity values can be found in Table S2. In parallel, the Kaiju tool was used to identify the different taxa present in the sample, which revealed a total of 2666 different taxa. Surprisingly, the vast majority of the identified taxa were of viral origin. As shown in Figures 2 and 3, 94.71% of the taxa significantly represented in this study were viral, highlighting the success of the combined enrichment procedure used for this study. A total of 4.6% of the taxa were found to be of bacterial origin and less than a total of 1% corresponded to cellular organisms (mostly fungi) and archaea. The full complement of taxa can be found in Table S3. The literature was scanned for studies that attempted to look at bat guano viromes applying a variety of viral enrichment methods [10,[28][29][30][31].    1, Wuhan Millipede virus 4, Helicoverpa armigera iflavirus, and Wenzhou tombu rus 11. The six most represented bacterial taxa included: Pseudomonas spp., Ac baumanii, Pseudomonas syringae, Burkholderia spp., and two unidentified Gamma teria. In the two underrepresented categories, one taxon corresponded to a fung and another one to a eukaryotic parasite (Plasmodium ovale wallikeri) that was p identified in African bats. Figure 3. Origin of the taxa with more than 0.1% representation in the data set. Each circle one of the taxa that represents >0.1% of the reads in the population. The color attribut circle identifies the domain of life to which the taxa belongs, as shown in the legend: pink yotes, blue for proteobacteria, terracotta for viruses, and light green for cellular organis The figure was generated using the Kaiju web server tool. Origin of the taxa with more than 0.1% representation in the data set. Each circle represents one of the taxa that represents >0.1% of the reads in the population. The color attributed to each circle identifies the domain of life to which the taxa belongs, as shown in the legend: pink for eukaryotes, blue for proteobacteria, terracotta for viruses, and light green for cellular organisms (fungi). The figure was generated using the Kaiju web server tool.
Determination of the proportion of total viral sequences allowed us to determine that our strategy enabled us to secure the highest number of viral reads for bat guano virome study to date, to our knowledge (Table 1). In an effort to try to define the representation of the different taxa present in the sample, we looked at the origin of the taxa that represented more than 0.1% of the total sample.  Figure 3 shows that the vast majority (66/74) of the most represented taxa were of viral origin, further supporting the success of the enrichment strategy used for this study. The ten most represented taxa corresponded to insect viruses. This is not surprising considering M. chiloensis' insectivorous diet. These taxa corresponded to Hubei diptera virus 21 two underrepresented categories, one taxon corresponded to a fungal species and another one to a eukaryotic parasite (Plasmodium ovale wallikeri) that was previously identified in African bats.
Considering the characteristics of the samples that were processed, including a pancoronavirus and panparamyxo virus PCR positive test and the proximity to humans, we interrogated whether any viruses in this family could be identified among the different taxa. Surprisingly, it was not possible to identify contings or taxa related to paramyxoviruses. However, both strategies enabled us to identify contigs with identities to alphacoronaviruses. On the one hand, the de novo assembly strategy allowed for identification of three contigs with high coverage that were found to have significant identity to bat coronaviruses. These contigs were used for mapping and reassembly of a full genome sequence that was then used for phylogenetic analyses (Figure 4). On the other hand, taxa classification using the Kaiju tool also identified a minimum of 12 contigs with significant identities to coronaviruses (Table S3). Out of these nine showed identities to alphacoronaviruses and five were likely to be of bat origin. These contigs were constructed with a minimum of 200 reads.
The three contigs generated by the de novo strategy were then used for mapping and reassembly in order to see if we would be capable of generating a full genome sequence. The full sequence was 28,323 nucleotides in length, suggesting the full genome sequence was generated. Phylogeny tree was built using 56 related sequences and the newly identified genome (Table S1). As seen in Figure 4, the newly identified genome sequence is grouped with six different bat alphacoronaviruses and a more distantly related Porcine epidemic diarrhea virus, further confirming the presence of this pathogen in this sample. Interestingly, most of the alphacoronaviruses of bat origin that were found in the same clade as the one identified in this study originated from China and one from Kenya and were isolated in insectivorous bats.

Discussion
Recently, several major outbreaks of emerging viral diseases, including the Hendra, Nipah, Marburg, and Ebola viruses, severe acute respiratory syndrome (SARS), and Middle East respiratory syndrome (MERS), as well as the current COVID-19 coronavirus pandemic, have been reported with increasing frequency [32,33]. Importantly, all these outbreaks have been linked to a suspected zoonotic origin from bat-borne viruses (or at least with bats playing a role as an intermediate host) [11]. Bats, the only mammal capable of powered flight display features that are unique among this animal class, a long lifespan relative to body size, a low rate of tumorigenesis, and an exceptional ability to carry diverse viruses without developing signs of clinical disease. Even though human-to-human transmission of SARS-CoV and MERS-CoV seems to be the most common type of event, some human cases recorded were occasionally the results of zoonotic transmission. This raised interest in the identification of animal reservoirs. In addition, an improved understanding of the natural history of these epidemics, including tracing viral evolution across different hosts, will allow us to improve how we tackle such public health challenges. Evidence suggests that almost all human coronaviruses have zoonotic origins or otherwise have a close relative that circulates in wild animals (bats) and domestic animals (camels and cattle). Besides SARS-CoVs and MERS-CoVs, numerous other CoVs have been detected in bats in Africa, Asia, Europe, and America and are classified into the genera alphacoronavirus and betacoronavirus. The metagenomic view on viruses (virome) unlocks the simultaneous detection of viruses in individual or pooled samples within a very short time. However, for samples with a high host material background, it is becoming apparent that purification procedures are necessary for NGS approaches to significantly increase the number of obtained viral sequences from biological samples.
In our study, we describe for the first time the virome of a previously metagenomically uncharacterized native bat species from Latin America: M. chiloensis. Importantly, the characterized colony lived in close proximity to humans, i.e., a church in the O'Higgins region. A combined approach for viral enrichment through benchtop centrifugations and efficient RNA depletion in an automated platform, the Magelia, enabled us to provide insights into the M. chiloensis virome with unprecedented resolution. To our knowledge, our strategy provided for the highest proportion of viral sequences in host "contaminated" samples Table 1.
Several strategies have been applied historically to enrich the viral content of guano samples used for virome studies. Centrifugation and filtration are the most commonly utilized strategies. Most studies incorporate these strategies to remove host cells or other micro-organisms. A nuclease-treatment step, where DNAse or RNAse will destroy exogenous nucleic acids but is not thought to affect nucleic acid protected by the viral capsid or envelope can be applied in addition. Finally, ultracentrifugation is also featured as a common method for the concentration of viruses from samples, applied in addition. Finally, ultracentrifugation is also featured as a common method for the concentration of viruses from samples. When looking at studies in the literature that attempted to enrich for similar samples as ours, we determined that our strategy achieved the highest level of enrichment when estimating the proportion of viral sequences. Studies that employed centrifugation, filtration, and Nuclease treatments were below our level of enrichment [28,29]. In Paskey's 2020 study, none of the aforementioned strategies were employed (or at least cited in the Materials and Methods), which could explain that this work reported the lowest level of enrichment. Interestingly, Salmier's study of 2017 incorporated the use of a benchtop instrument for automated nucleic acid extraction, the NucliSENS easyMAG (Biomérieux). The instrument's proprietary magnetic silica particle manipulation technology did not seem to impact viral nucleic acid enrichment, which is not surprising as it is not specifically targeted for this. Nucleic acid extraction methods utilized for these studies were similar except for the latter study. This suggests that this didn't impact the level of enrichment. Our strategy included centrifugation, filtration, nuclease treatment, and in addition to this, negative enrichment through probe hybridization to select for host and microbial nucleic acids in an automated benchtop instrument. It is likely that Paskey also used a similar approach (albeit, without automation), as he used a previous version of the Illumina kit we used that likely includes probe negative enrichment. This is not specified in Table 1 as it is not explicit in the publication's material and methods section. Despite this, it achieved lesser viral enrichment. The key to our approach resides in miniaturization of the reaction which improves kinetics paired with highly precise manipulation of magnetic beads, improving all cleanup steps. Our study shows that enrichment of such complex samples for virome studies can be improved utilizing a combination of traditional techniques including centrifugation, filtration, and nuclease treatment paired with molecular negative enrichment utilizing a mixture of probes. Automation, miniaturization, and precise bead handling thanks to a patented magnetic tweezer further improve enrichment and impact the resolution of the obtained data.
On the first level of resolution, the data generated in this study paints a detailed picture of the M. chiloensis virome. Unsurprisingly, most identified viruses are known to be associated with diverse insect species. This will ultimately enable us to gain insights into the diet of this bat species, which will in turn improve our knowledge on the role of M. chiloensis in the ecosystem and clarify its role as a pest-controlling mammal. Additionally, the resolution of the data obtained here could help indirectly survey insect population health, for example by identifying viruses affecting pollinator species. Although the dataset generated here is sufficient for this endeavor, this is not the main focus of this study and these analyses, although ongoing, are outside of the scope of this manuscript.
On a second level, the identification of diverse pathogenic viruses of zoonotic interest was made possible by the degree of enrichment of our study. Due to the fact that the sample tested positive for a pancoronavirus PCR test, this family of viruses was of particular interest. De novo assembly in parallel with taxonomic classification allowed for the identification of a novel alphacoronavirus sequence closely related to PEDV. A bat coronavirus sequence from Tadarida brasiliensis was linked with PEDV in Brazil [34], which reinforces the plausibility of our results.
In China, an alphacoronavirus, bat HKU2-like genetically distinct but clinically similar to PEDV was identified in diarrheic pigs and was designated swine acute diarrhea syndrome coronavirus (SADS-CoV) [35][36][37][38]; however, this genome sequence was phylogenetically far from alphacoronavirus detected in M. chiloensis. Furthermore, human alphacoronavirus 229E and Nl3 were not related to the alphacoronavirus genome identified in this study. This, however, should not minimize the relevance of determining its potential to jump to humans. Recently, a coronavirus isolated in Indonesian bats showed the capability of infecting human cells [39], highlighting once again the importance of sustaining this type of epidemiological surveillance. In addition to this, it should be noted that we did not detect sequences associated with any Rabies virus or coronavirus from the beta genera, underlining the absence of known zoonotic threats in this population. Our data improves the understanding of the viruses currently found in native chilean bats, highlighting the need to test the zoonotic potential of the virus identified through this study.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/v14020202/s1, Table S1: Accession numbers for the viral sequences used for alignment/phylogenetic analyses in this study; Table S2: Significant identity values for generated contigs blasted in NCBI's nr database; Table S3


Association Between GSDMB Gene Polymorphism and Cervical Cancer in the Northeast Chinese Han Population

Objective: The purpose of this study was to investigate the relationship between GSDMB gene polymorphism and genetic susceptibility to cervical cancer in the Han population in Northeast China. Methods: In this case–control study, the genotypes and alleles of rs8067378 in the GSDMB gene were analyzed by multiplex polymerase chain reaction (PCR) and next-generation sequencing methods in 482 cervical cancer (CC) patients, 775 cervical squamous intraepithelial lesion (SIL) patients, and 495 healthy women. The potential relationships between the SNP of the GSDMB gene with SIL and CC were analyzed by multivariate logistic regression analysis combined with 10,000 permutation tests. Results: In the comparison between the SIL group and the control group, the genotype and allele distribution frequencies of rs8067378 SNP of the GSDMB gene were statistically significant (p = 0.0493 and p = 0.0202, respectively). The allele distribution frequencies of rs8067378 were also statistically significant in the comparison between high-grade cervical squamous intraepithelial lesion (HSIL) and low-grade cervical squamous intraepithelial lesion (LSIL) groups with control group ( p = 0.0483 and p = 0.0330, respectively). Logistic regression analysis showed that after adjusting for age, the rs8067378 SNP of the GSDMB gene was significantly associated with the reduced risk of SIL under the dominant model (p = 0.0213, OR = 0.764, CI = 0.607–0.961) and the additive model (p = 0.0199, OR = 0.814, and CI = 0.684–0.968), and its mutant gene G may play a role in the progression of healthy people to LSIL and even HSIL as a protective factor. However, there was no significant association between cervical cancer and its subtypes with the control group (p > 0.05). After 10,000 permutations, there was still no correlation that has provided evidence for the accuracy of our study. Conclusion: The results of this study showed that rs8067378 single nucleotide polymorphism of the GSDMB gene may reduce the risk of SIL and protect the susceptibility to cervical precancerous lesions in the Northeast Chinese Han population, but it has no significant correlation with the progression of cervical cancer.

INTRODUCTION
Cervical cancer is the fourth most frequently diagnosed cancer and the fourth leading cause of cancer death in women. It is estimated that by 2020, there were 604,000 new cases and 342,000 deaths in the world (Sung et al., 2021). The progression of cervical cancer is a complex process with multiple steps and factors. In addition to the risk factors such as long-term chronic inflammation caused by high-risk human papillomavirus (HPV) persistent infection, longterm oral contraceptives, and fertility and sexual life disorder, the occurrence of cervical cancer also has obvious individual genetic susceptibility (Nkfusai et al., 2019;Olusola et al., 2019;Sadri Nahand et al., 2020;van der Waal et al., 2020;Yang et al., 2020). Single nucleotide polymorphism (SNP) is the most common form of genetic variation, which widely exists in the human genome and is related to the susceptibility to a variety of cancers such as cervical cancer. Therefore, the study of SNP will also contribute to the early diagnosis and treatment of cervical cancer (Shastry, 2009;Tan, 2017).
Gasdermin B (GSDMB) is a member of the GSDM protein superfamily that exists only in vertebrates. The GSDMB gene, located on chromosome 17q21.2, has a conserved gasdermin domain, a potential nuclear localization signal, and two nuclear receptor-binding motifs LXX(L/I)L. It is 1518 bp length and encodes 411 amino acids. It is often expressed in normal esophageal and gastrointestinal epithelium and bronchial epithelium of asthmatic lungs (Sun et al., 2008;Zheng et al., 2020). GSDMB has six splice variants, and GSDMB-1 is the longest and most important subtype which is expressed in human cancer cell lines (Feng et al., 2018).
Recently, several studies have identified that the expression of GSDMB has a potential relationship with the occurrence and development of tumors. GSDMB can regulate its lipid binding and pore-forming activity through different mechanisms of the intramolecular domain and then participate in pyrocytosis (Chen et al., 2019;Li et al., 2020). Pyrocytosis is a proinflammatory form that regulates cell death. Its essence is a    CC, cervical cancer; SIL, cervical squamous intraepithelial lesion; OR, odds ratio; CI, confidence interval. P: p-value (CC vs. SIL ); p-values <0.05 are shown in bold. P': age was adjusted. P": the p-value was calculated using 10,000 permutations for each model to correct the multiple tests. cascade amplification of inflammatory response. It often participates in the body's immune defense and plays a "double-edged sword" role in the occurrence and development of tumor (Li et al., 2020). In addition, the level of sulfatide on the apical surface of many tumor epithelial cells is significantly increased, and the metastatic ability of tumor cells is related to the specific binding of GSDMB to sulfatide (Chao et al., 2017). Increasingly, a large number of case-control studies have also suggested that the GSDMB gene polymorphism is associated with susceptibility to many autoimmune diseases and malignant tumors such as cervical cancer (Hergueta-Redondo et al., 2016;Wiemels et al., 2018;Molina-Crespo et al., 2019;Li et al., 2020;Li et al., 2021). At present, the relationship between GSDMB gene polymorphism and the risk of cervical cancer in the Han population in Northeast China has not been reported. Therefore, this study aimed to explore the relationship between rs8067378 polymorphism of the GSDMB gene and the risk of cervical cancer in the Northeast Chinese Han population. Since cervical squamous intraepithelial lesions, as a form of precancerous lesions, reflect the continuous development of cervical cancer (de Rycke et al., 2020), we will also study whether rs8067378 polymorphism is associated with cervical squamous intraepithelial lesions in Han women in Northeast China. HSIL, high-grade cervical squamous intraepithelial lesion; LSIL, low-grade cervical squamous intraepithelial lesion; OR, odds ratio; CI, confidence interval. P: p-value (HSIL, vs. LSIL ); p-values <0.05 are shown in bold. P': age was adjusted. P": the p-value was calculated using 10,000 permutations for each model to correct the multiple tests.

Control Group
The control group consisted of 496 healthy people from the physical examination center of the same hospital in the same period. Inclusion criteria: there was no abnormality in the ThinPrep cytologic test (TCT), no cancer, or family history of cancer. Exclusion criteria: all gynecological diseases and surgical history of gynecological diseases, hypo-immunity and immune diseases, cervical surgery, skin or genital condyloma, and other cancer history, etc. All subjects signed informed consent, and all studies in this report were approved by the Ethics Committee of the Second Affiliated Hospital of Harbin Medical University.

Candidate SNP Selection
Based on the data from the previous literature and combined with the characteristics of the East Asian population in the dbSNP database, the SNP rs8067378 (A > G) at 17q12 identified by previous cervical cancer GWAS (Shi et al., 2013) was selected for the current study.

Extraction and Genotyping of DNA
After fasting for 12 h, 1 ml peripheral venous blood samples were collected and put into a 2% EDTA-Na₂ anticoagulant tube. The samples were numbered and stored in an −80°C refrigerator until DNA extraction. The genomic DNA of all subjects was extracted from peripheral venous blood samples for genotyping according to the standard steps of the TIANamp Genomic DNA Kit (Tiangen Biotech, Beijing, China). All genotyping experiments for the selected SNP (rs8067378) were performed using two rounds of multiplex PCR combined with next-generation sequencing methods by the Shanghai Biowing Applied Biotechnology Company (Chen et al., 2016). Primer3 online software (v0.4.0) was used to amplify the primer sequences. Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 860727 5 PCR was performed for SNP analysis with the following primers: rs8067378 (sense, 5′-GTTGACAGTCAAAACAAAAACCTG-3′). After two rounds of multiplex PCR, the PCR products were mixed into the centrifugal tube, then sealed with the sealing membrane, and placed overnight. The mixture was purified using a TIANgel Midi Purification Kit (Tiangen Biotech, China). The purified PCR products were then pairedend sequenced (2 × 150 bp) by the Illumina HiSeq XTen platform according to the manufacturer's instructions. The read data were compared to the human reference genome using the Burrows-Wheeler (BWA, v0.7.12) (Slater et al., 2009), and SAMtools (v0.1.19) (Li, 2011) was used for SNP calling and genotyping. In total, ninety-one samples were randomly selected for blind DNA replication to control the quality of genotyping.

Physiological and Biochemical Index Collection and Analysis
The general information and clinical data of each subject in the case group were recorded in detail: age, menstrual history, smoking and drinking history, past disease history, operation history, family history of tumor, relevant laboratory examination results, histopathological diagnosis, clinically confirmed diagnosis, and relevant gynecological examination (ultrasound or colposcopy), etc.
In addition, the general situation and some clinical data of the healthy control population in the physical examination center during the same period were recorded: age, ThinPrep cytologic test (TCT), and HPV, etc. Moreover, due to the influence of traditional Chinese culture, the subjects in this study almost refused to answer questions about sexual life history, so this part of the data could not be collected.
Squamous cell carcinoma (SCC) antigen, carcino-embryonic antigen (CEA), alpha-fetoprotein (AFP), CA125, D-dimer, and other biochemical indicators were detected by the laboratory of the Second Affiliated Hospital of Harbin Medical University. The histopathological diagnosis of cervical tissue samples was confirmed by the pathology department of the Second Affiliated Hospital of Harbin Medical University.

Allele, Genotypes, and Statistics
SPSS 21.0 software was used to analyze the clinical characteristics of the case group and the control group. Student's t-test was used to compare the continuous variables among the groups. The Pearson chi-square test was used to compare the categorical variables between groups (the Fisher exact test was used when the expected value was less than 5). Plink1.9 single nucleotide polymorphism software was used for genotyping analysis. The Hardy-Weinberg equilibrium was analyzed. We established the following three models: dominant model: major allele homozygous vs. heterozygous + minor allele homozygous; Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 860727 6 recessive model: major allele homozygous + heterozygous vs. minor allele homozygous; and additive model: major allele homozygous vs. heterozygous vs. minor allele homozygous. The 95% confidence interval (CI) and odds ratio (OR) of each genetic model were calculated. Multivariate logistic regression analysis was performed to adjust age to evaluate the association of SNP with cervical cancer, cervical intraepithelial neoplasia, and their subtypes. Moreover, the chi-square test combined with logistic regression analysis was used to analyze the relationship between genotypes and alleles with clinical parameters. All statistical analyses showed that p < 0.05 had statistical significance.

SNP Genotype and Quality Controls
The SNP of the GSDMB gene was rs8067378, and its genotype distributions in the CC and SIL groups conform to the Hardy-Weinberg equilibrium (p = 0.9108 and p = 0.9508, respectively). In addition, we controlled the quality of 73 samples in this experiment, and the accuracy rate is 98.5%, and the error rate is within a reasonable range, which ensures the reliability and repeatability of the follow-up research results.

Clinical Characteristics of the Study Population
The clinical characteristics of all CC patients, SIL patients, and healthy controls in this study were statistically analyzed, as shown in Table 1. The results showed that compared with the healthy control group, the average age of CC and SIL patients was old, and the comparison between the control group and CC and SIL groups was statistically significant (p < 0.0001). The analysis of other clinical characteristics of the CC group and SIL group showed that there were significant differences in menarche age and amenorrhea between the two groups (p < 0.05). In addition, smoking was significantly correlated between the two groups (p = 0.0091) while drinking was not (p = 0.4952).

Analysis of Genotypes and Allele Distribution Frequencies
The genotypes and allele distribution frequencies of GSDMB gene polymorphism (rs8067378) in the CC group and its subgroups with the healthy control group are shown in Table 2. The results of the chi-square test showed that the genotypes and allele frequencies of rs8067378 single nucleotide polymorphism of the GSDMB gene had no significant difference between the CC group and healthy control group (p > 0.05), and there was no significant difference between the SCC and AUC groups with healthy control group (p > 0.05). Logistic regression analysis showed that after adjusting for age, rs8067378 single nucleotide polymorphism was not significantly associated with cervical cancer and its subtypes (p > 0.05). After 10,000 permutations, the results were still not statistically significant. The consistency of 10,000 permutations further supports the reliability of our results. Table 3 shows the comparison of genotypes and allele distribution frequencies between the SIL group and its subgroups with healthy control group. The results of the chi-square test showed that the genotype frequencies of rs8067378 single nucleotide polymorphism of the GSDMB gene were statistically significant between the SIL group and healthy control group (p = 0.0493), and the A/G allele frequency was also statistically significant between the two groups (p = 0.0202). Compared with the healthy control group (28.6%), the G mutation of rs8067378 accounted for 32.5% in the HSIL group and 34.9% in the LSIL group, and the A/G allele frequency was significantly correlated between the HSIL group and LSIL group with the healthy control group (p = 0.0483 and 0.0330, respectively). Logistic regression analysis showed that after adjusting for age, rs8067378 SNP was statistically significant between the SIL group and its subgroup (HSIL group) with the healthy control group in the dominant model and additive model (p = 0.0213 and 0.0356 in the dominant model, and p = 0.0199 and 0.0406 in the additive model, respectively). The difference between the LSIL group and healthy control group was statistically significant in the additive model (p = 0.0308). It can be seen that under the additive model, rs8067378 plays a protective role in the progression of healthy people to LSIL and even HSIL (P2 = 0.0308, OR = 0.740, CI = 0.564-0.973; P1 = 0.0406, OR = 0.827, CI = 0.689-0.992). This correlation still exists after 10,000 permutations. Table 4 shows the comparison of genotypes and allele distribution frequencies of rs8067378 SNP of the GSDMB gene between the CC group and SIL group. The results show that there is no significant difference in genotypes and allele distribution frequencies of rs8067378 SNP between the CC group and SIL group (p > 0.05). Logistic regression analysis showed that there was no significant difference in rs8067378 SNP between the CC group and SIL group after adjusting for age (p > 0.05). After 10,000 permutations, the results were still not statistically significant.
Finally, we made intra-subgroup comparisons. Tables 5, 6 summarized the comparisons of genotypes and allele distribution frequencies of rs8067378 SNP of the GSDMB gene between the SCC group and the AUC group, and the HSIL group and the LSIL group, respectively. The results show that there is no significant difference in genotypes and allele distribution frequencies of rs8067378 SNP between the two intra-subgroups (p > 0.05). After adjusting for age, there is still no significant difference among the genetic models (p > 0.05). After 10,000 permutations, the results were still not statistically significant, which has provided evidence for the accuracy of our study. In addition, in order to more clearly and definitely show the comparison of genotypes and allele distribution frequencies between groups, we drew three graphs: Figure 1 shows the results of chi-square tests, Figure 2 shows the results of logistic regression analysis after adjusting for age, and Figure 3 shows the results after 10,000 permutation tests.

Analysis of Genotypes and Allele Frequencies With Clinical Parameters
In the CC group and SIL group, the comparisons of genotypes and allele distribution frequencies of GSDMB gene rs806778 polymorphism with clinical parameters are shown in Tables 7-10. It is worth noting that in patients with cervical cancer, allele frequencies are significantly correlated with amenorrhea (p = 0.0328). After logistic regression analysis, the correlations still exist (p = 0.0330). In addition, no significant correlation was observed between genotypes and allele frequencies with other clinical parameters (p > 0.05). Moreover, there was also no significant correlation between genotypes and allele frequencies with clinical parameters in patients with cervical squamous intraepithelial lesions (p > 0.05).

DISCUSSION
Cervical cancer is characterized by a high incidence rate and high mortality rate, and it is a major public problem affecting the health of middle-aged women, especially in countries with Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 860727 8 low resources (Hu et al., 2018;Arbyn et al., 2020;Zhang et al., 2020). Persistent infection of high-risk HPV is the most important risk factor associated with cervical cancer (Faridi et al., 2011;Paaso et al., 2019). However, studies have shown that other risk factors associated with HPV may also play important roles in the pathogenesis of cervical cancer, such as immune or genetic factors. SNP has become an important biomarker for locating cancer, such as cervical cancer (Shastry, 2009). Therefore, this case-control study was conducted to research the association between GSDMB gene polymorphism and the risk of cervical cancer in the Han population in Northeast China. In this study, we investigated the relationship between rs8067378 SNP of the GSDMB gene and susceptibility to cervical cancer in Han population in Northeast China. According to the histopathological classification, CC patients were further divided into the SCC group and AUC group, and the SIL group was divided into the HSIL group and LSIL group. Then, intragroup and intergroup comparisons were performed.
Our results replicated previous studies and showed that the average age of the CC group was 49.33 ± 9.29 years, and the average age of the SIL group was 42.22 ± 9.68 years. This is consistent with the conclusion that the peak of CC in the Chinese mainland is 40-60 years old and the peak of SIL is 30-50 years old (Wang et al., 2018;Liu et al., 2020). This may be related to the fact that women in the 30-40 age group are in the period of sexual activity and childbearing, while most of the women after the age of 50 are postmenopausal (Misra et al., 2018).
In this case-control study, we investigated the relationship between rs8067378 polymorphism of the GSDMB gene and the risk of CC in the Han population in Northeast China. The results showed that there was no significant correlation between rs8067378 polymorphism with CC and its subtypes. These are different from the existing research studies. A genome-wide association study on cervical cancer in the Chinese Han population proved that rs8067378 SNP is a susceptible site for cervical cancer (Shi et al., 2013). Miura et al. (2016) showed for the first time that there was a significant correlation between rs8067378 SNP and invasive cervical cancer and further proved that Japanese women with the GG genotype were at high risk of invasive cervical cancer. Lutkowska et al. (2017) studied the relationship between rs8067378 SNP and stage III and IV of cervical cancer in the Polish population, in which the G allele plays a role in the diffusion of tumor cells to adjacent tissues, indicating that rs8067378 SNP increases the risk of occurrence and development of cervical cancer. Such findings have also been confirmed in some parts of southern China. Based on previous studies, the specific mechanism of GSDMB in the occurrence and development of cervical cancer is not clear. Some studies suppose that GSDMB may participate in the regulation of the estrogen-estrogen receptor-target gene expression pathway under the action of endogenous and exogenous estrogen (Sun et al., 2008). Studies have shown that GSDMB containing two nuclear receptor-binding motifs can be used as a nuclear receptor co-activator, recruited in the estrogen receptor to form a complex, provide related enzymatic activity and scaffolding function, promote the high expression of E6/E7 oncogene in patients with high-risk HPV persistent infection, and cause uncontrolled cell proliferation, cervical intraepithelial lesions, and even cervical cancer (Sun et al., 2008). In addition, these differences may be related to the sample size and case determination protocol. Also, the genetic background and lifestyle may also change the genetic susceptibility of the GSDMB gene, which explain the differences in the intensity of research correlation between different regions and races. Therefore, we should repeat the study with a larger sample, to determine the relationship between GSDMB gene polymorphism with SCC and AUC and even other special types of CC.
GSDMB can participate in pyrocytosis and the production of inflammatory factors, and studies found that GSDMB can be highly expressed in cervical inflammatory lesions and precancerous lesions (Sun et al., 2008). Therefore, in order to further study the role of GSDMB gene rs8067378 polymorphism in the development of cervical cancer in healthy women, we also analyzed the susceptibility of SIL. Surprisingly, under the dominant and additive models, the rs8067378 polymorphism of the GSDMB gene is significantly correlated with SIL and its subtypes, and rs8067378 single nucleotide polymorphism may reduce the risk of SIL and protect the susceptibility to SIL. This may be related to the dual role of GSDMB-mediated pyrocytosis in tumor, and the occurrence of pyrocytosis is concentrationdependent on GSDMB. Studies have pointed out that (Kayagaki et al., 2015;Tsuchiya, 2020) macrophages and monocytes lacking GSDMB almost completely lose the transmission ability of lipopolysaccharide (LPS) in the cytoplasm and will not activate inflammatory caspases, so they will not induce pyrocytosis. In the process of pyrocytosis, caspases are the inducer and GSDMB is the substrate. Precursor caspases are activated by a variety of inflammatory bodies such as the NODlike receptor (NLR) and are absent in melanoma 2 (AIM2). Active caspases bind to the cutting site of GSDMB. This process releases the N-terminal and C-terminal, which are inserted into the cell membrane to cause pyrocytosis and release a large number of inflammatory factors such as tumor necrosis factor (TNF), interleukin-1β (IL-1β), and interferon-α (IFN-α). This process can not only directly kill pathological and tumor cells but also further recruit immune cells, so as to expand the inflammatory response. Therefore, some scholars believe that active caspases cutting GSDMB-induced pyrocytosis is a positive feedback mechanism in the immune process, so that most newborn tumor cells are eliminated by the process of pyrocytosis before the formation of tumor. It plays a role in preventing the further development of precancerous lesions (Dinarello, 2009;Guo et al., 2015;Xia et al., 2019;Zhou and Fang, 2019;Ruan et al., 2020). This is consistent with the research results of Zhang WH. The study has shown that after CC cells are infected by HPV, AIM2 can play a tumor inhibitory role by inducing pyrocytosis (So et al., 2018).
The study confirmed for the first time that GSDMB gene polymorphism was significantly correlated with cervical squamous intraepithelial lesions, but not with cervical cancer in the Han population in Northeast China. The studies of Miura et al. and Lutkowska et al. have confirmed that GSDMB gene polymorphism was significantly associated with cervical cancer. Compared with previous existing studies, the sample size in our study was relatively large, including 482 cases of cervical cancer, 775 cases of cervical squamous intraepithelial lesions, and 495 healthy controls. Moreover, our study was analyzed from multiple perspectives such as cervical cancer subgroups and cervical precancerous lesions. A large sample size is very important for exploring and verifying the potential mechanism of cervical cancer and precancerous lesions in Chinese women and will also play an important role in the research of different pathological types of cervical cancer and precancerous lesions in the future. In addition, we also added the comparison of genotypes and alleles with clinical parameters to understand more clinically related risk factors. However, this study also has some limitations. First, this study adopts a case-control


Metabolic Changes in Larvae of Predator Chrysopa sinica Fed on Azadirachtin-Treated Plutella xylostella Larvae

Biological control is a key component of integrated pest management (IPM). To suppress pests in a certain threshold, chemical control is used in combination with biological and other control methods. An essential premise for using pesticides in IPM is to ascertain their compatibility with beneficial insects. Chrysopa sinica (Neuroptera: Chrysopidae) is an important predator of various pests and used for pest management. This study was intended to analyze metabolic changes in C. sinica larvae after feeding on azadirachtin-treated Plutella xylostella (Lepidoptera, Plutellidae) larvae through a non-targeted LC–MS (Liquid chromatography–mass spectrometry) based metabolomics analysis. Results showed that C. sinica larvae did not die after consuming P. xylostella larvae treated with azadirachtin. However, their pupation and eclosion were adversely affected, resulting in an impairment in the completion of their life cycle. Feeding C. sinica larvae with azadirachtin-treated P. xylostella larvae affected over 10,000 metabolites across more than 20 pathways, including the metabolism of amino acids, carbohydrates, lipid, cofactors, and vitamins in C. sinica larvae, of which changes in amnio acid metabolism were particularly pronounced. A working model was proposed to illustrate differential changes in 20 metabolites related to some amino acid metabolisms. Among them, 15 were markedly reduced and only five were elevated. Our results suggest that azadirachtin application may not be exclusively compatible with the use of the predator C. sinica for control of P. xylostella. It is recommended that the compatibility should be evaluated not only based on the survival of the predatory insects but also by the metabolic changes and the resultant detrimental effects on their development.

Introduction
Integrated pest management (IPM) is a coordinated process using multiple methods, such as biological, chemical, cultural, mechanical, physical, and pest resistant or tolerant varieties for optimizing control of pests in an ecologically and economically sound manner [1]. Among them, chemical control and biological control are two common methods used either separately or in combination for managing insect pests during crop production. However, the use of pesticides may have adverse effects on non-target organisms, including predatory insects [2]. In order to minimize these problems and maintain sustainable control of insect pests, botanical pesticides are considered to be attractive alternatives for pest management. Compared with traditional synthetic insecticides, botanical pesticides have better eco-toxicological properties, including low toxicity, rapid degradation, and little impact on the environment, which makes them a suitable choice for pest control [3,4].

Azadirachtin Activities against C. sinica
The larvae of C. sinica fed with P. xylostella larvae that consumed azadirachtin-treated leaves (T) showed no significant growth differences compared to those fed P. xylostella control (CK) as there were no larval mortalities between T and CK before pupation (data not shown). After stopping preying, mature larvae began to pupate. The head and tail of the larvae gradually curled together, and the tail drew silk to make cocoons (Figure 1(c-4)). However, there were 18% of the larvae in the T treatment that could not curl up and could not draw silk to make a cocoon (Figure 1(c-1-c-3)). When touching the larvae with a small brush, they still twisted but died in a few days. After 20 to 30 days, the C. sinica adults emerged from the pupae. In the T treatment, 24% of adults were deformed as they were unable to extend their wings and/or had malformed abdomen (Figure 1(d-1)), which were regarded as the failure to eclosion. Adults with fully extended wings and no growth defects were considered as successful eclosion (Figure 1(d-2,d-3)). Thus, the pupation and eclosion of C. sinica larvae were significantly affected by T treatment. As shown in Figure 2, the proportion of larvae undergoing pupation from T treatment was significantly lower at 82.00 ± 3.06% comparted to 100.00% in the CK (p < 0.01), and their eclosion was 76.00 ± 1.15% in the T treatment against 98.00 ± 1.48% of CK (p < 0.001).   The process of monitoring P. xylostella larvae feeding on cabbage leaves and C. sinica larvae ingesting azadirachtin-treated P. xylostella larvae and subsequently their pupation and eclosion. A larva feeding on a cabbage leaf treated with azadirachtin (a). A C. sinica larva ingesting a P. xylostella larva (b). C. sinica larvae underwent pupation from curling to the formation of cocoon (c-1-c-4). The emergence of C. sinica adult: a deformed adult (d-1), normal adult (d-2), and normal adult after removing wings to show normal abdomen (d-3).

Figure 2.
The proportion of mature C. sinica larvae underwent pupation (a) and eclosion (b) after ingesting azadirachtin-treated P. xylostella larvae. Data were expressed as the mean ± S.E. and ** and *** indicate significant differences at p < 0.01 and p < 0.001 levels based on Tukey's HSD test.

Metabolic Profiles Analyzed by LC-MS
The unsupervised PCA was used to check the quality of the data from the LC-MS analyses. In ESI+ mode, the PC1 and PC2 explained 30.9% and 11.4% of the total variance of all samples. In ESI− mode, the PC1 and PC2 explained 33.9% and 10.3% of the total variance. The supervised PLS-DA was performed to identify the metabolites responsible for the separation between CK and azadirachtin treatments (T). Results showed that in the ESI+ mode, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the PLS-DA model were 0.59, 0.753, 0.369, and 0.269 (Figure 3a), respectively; in the ESI− mode, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the PLS-DA model were 0.728, 0.74, 0.536, and 0.276 (Figure 3b), respectively. Based on the OPLS-DA model, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the ESI+ mode were 0.59, 0.753, 0.337, and 0.269 (Figure 4a), respectively. In the ESI− mode, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the ESI− mode were 0.728, 0.74, 0.512 and 0.276 (Figure 4b), respectively. ingesting azadirachtin-treated P. xylostella larvae. Data were expressed as the mean ± S.E. and ** and *** indicate significant differences at p < 0.01 and p < 0.001 levels based on Tukey's HSD test.

Metabolic Profiles Analyzed by LC-MS
The unsupervised PCA was used to check the quality of the data from the LC-MS analyses. In ESI+ mode, the PC1 and PC2 explained 30.9% and 11.4% of the total variance of all samples. In ESI− mode, the PC1 and PC2 explained 33.9% and 10.3% of the total variance. The supervised PLS-DA was performed to identify the metabolites responsible for the separation between CK and azadirachtin treatments (T). Results showed that in the ESI+ mode, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the PLS-DA model were 0.59, 0.753, 0.369, and 0.269 (Figure 3a), respectively; in the ESI− mode, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the PLS-DA model were 0.728, 0.74, 0.536, and 0.276 (Figure 3b), respectively. Based on the OPLS-DA model, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the ESI+ mode were 0.59, 0.753, 0.337, and 0.269 (Figure 4a), respectively. In the ESI− mode, the R 2 X, R 2 Y, Q 2 Y, and RMSEE values in the ESI− mode were 0.728, 0.74, 0.512 and 0.276 (Figure 4b), respectively.

Changes in Metabolites and Metabolic Pathways of Differentially Abundant Metabolites
Representative LC-MS total ion chromatograms (TICs) of C. sinica larvae tissue samples are shown in Figure 5. The shape and quantity of peaks between the T and CK treatments varied greatly. Approximately 13,672 and 10,947 metabolite peaks were deconvoluted in ESI+ and ESI− modes of LC-MS, respectively. The ESI+ usually detects N, O, and S-containing species and also some specific hydrocarbons, such as isoprene, terpenes, and aromatics as protonated neutral MH+, whereas the ESI− detects acid including carboxylic acids RCOOH and inorganic acids and hydrosulfides as deprotonated neutral [M-H] − [42]. A total of 3210 and 2026 remaining peaks in ESI+ and ESI− modes in LC-MS were further annotated using references in existing databases, respectively. After the exogenous compounds in LC-MS were removed, the differentially abundant metabolites were selected according to the VIP values from the OPLS-DA model (VIP ≥ 1) and the corrected p values from t-test (p < 0.05). There were 778 compounds in the ESI+ model, of which 357 were

Changes in Metabolites and Metabolic Pathways of Differentially Abundant Metabolites
Representative LC-MS total ion chromatograms (TICs) of C. sinica larvae tissue samples are shown in Figure 5. The shape and quantity of peaks between the T and CK treatments varied greatly. Approximately 13,672 and 10,947 metabolite peaks were deconvoluted in ESI+ and ESI− modes of LC-MS, respectively. The ESI+ usually detects N, O, and S-containing species and also some specific hydrocarbons, such as isoprene, terpenes, and aromatics as protonated neutral MH+, whereas the ESI− detects acid including carboxylic acids RCOOH and inorganic acids and hydrosulfides as deprotonated neutral [M-H] − [42]. A total of 3210 and 2026 remaining peaks in ESI+ and ESI− modes in LC-MS were further annotated using references in existing databases, respectively. After the exogenous compounds in LC-MS were removed, the differentially abundant metabolites were selected according to the VIP values from the OPLS-DA model (VIP ≥ 1) and the corrected p values from t-test (p < 0.05). There were 778 compounds in the ESI+ model, of which 357 were upregulated and 421 were downregulated. In ESI− mode, there were 391 compounds: 180 were upregulated and 211 were downregulated.   Table 1 shows some related metabolic pathways with representative differentially abundant metabolites and their upregulation and downregulation. Amino acids, carbohydrates, bile, lipids, membrane transports, cofactors and vitamins were primary metabolites that were affected in C. sinica larvae after ingestion of azadirachtin-treated P. xylostella larvae. The enrichment of pathways is presented in Figure 6. Rich factor refers to the ratio of the numbers of differentially abundant metabolites annotated in this pathway to the numbers of all metabolites annotated to the same pathway. The greater the rich factor, the greater the pathway enrichment. The p value is another parameter for enrichment with a range from 0 to 1, the closer to 0, the more significance of the enrichment. As shown in Figure 6, biotin metabolism was significantly enriched, but the numbers of metabolites were much lower. Tryptophan metabolism was also significantly enriched with relatively higher metabolite numbers. Arginine and proline metabolism had a Rich factor of 0.33 with rather higher metabolite numbers. Lysine degradation had comparable Rich factors and also higher metabolite numbers. Additionally, beta-alanine metabolism was significantly relevant pathway of C. sinica larvae that were affected by ingestion of azadirachtin-treated P. xylostella larvae.  Table 1 shows some related metabolic pathways with representative differentially abundant metabolites and their upregulation and downregulation. Amino acids, carbohydrates, bile, lipids, membrane transports, cofactors and vitamins were primary metabolites that were affected in C. sinica larvae after ingestion of azadirachtin-treated P. xylostella larvae. The enrichment of pathways is presented in Figure 6. Rich factor refers to the ratio of the numbers of differentially abundant metabolites annotated in this pathway to the numbers of all metabolites annotated to the same pathway. The greater the rich factor, the greater the pathway enrichment. The p value is another parameter for enrichment with a range from 0 to 1, the closer to 0, the more significance of the enrichment. As shown in Figure 6, biotin metabolism was significantly enriched, but the numbers of metabolites were much lower. Tryptophan metabolism was also significantly enriched with relatively higher metabolite numbers. Arginine and proline metabolism had a Rich factor of 0.33 with rather higher metabolite numbers. Lysine degradation had comparable Rich factors and also higher metabolite numbers. Additionally, beta-alanine metabolism was significantly relevant pathway of C. sinica larvae that were affected by ingestion of azadirachtin-treated P. xylostella larvae. KEGG is the major public pathway-related database that includes not only genes but metabolites. † Pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways in differential metabolites comparing with the whole background. The calculating formula is as follows: Here N is the number of all metabolites that with KEGG annotation, n is the number of differential metabolites in N, M is the number of all metabolites annotated to specific pathways, and m is number of differential metabolites in M. ‡ A variable importance in projection score of OPLS model was applied to rank the metabolites that best distinguished between two groups. § Means mass-to-charge ratio. Metabolites 2022, 12, x FOR PEER REVIEW 8 of 20 Figure 6. Metabolome map of significant metabolic pathways in C. sinica larvae affected by the ingestion of azadirachtin-treated P. xylostella larvae (pathway enrichment). Rich factor refers to the ratio of the number of annotated to this pathway in the differential metabolites to the number of annotated to this pathway in all metabolites. A larger rich factor indicates a higher degree of enrichment. p values range from 0 to 1, the closer to 0, the more significance of the enrichment. Figure 6. Metabolome map of significant metabolic pathways in C. sinica larvae affected by the ingestion of azadirachtin-treated P. xylostella larvae (pathway enrichment). Rich factor refers to the ratio of the number of annotated to this pathway in the differential metabolites to the number of annotated to this pathway in all metabolites. A larger rich factor indicates a higher degree of enrichment. p values range from 0 to 1, the closer to 0, the more significance of the enrichment.

Metabolite Changes Adversely Affected C. sinica Development
The present study showed that no mortality occurred in C. sinica larvae after consuming azadirachtin-treated P. xylostella larvae. However, their pupation and eclosion were significantly affected. Compared with the control treatment, 18% and 24% of the mature larvae were unable to perform pupation and eclosion, respectively (Figures 1 and 2). These results suggest that C. sinica larvae were able to obtain needed nutrients through digesting P. xylostella larvae to sustain their growth and even tolerate ingested azadirachtin, but the azadirachtin adversely affect C. sinica metabolism. As shown in Table 1 and Figure 6, over 10,000 metabolites across more than 20 pathways, including the metabolism of carbohydrates, lipid amino acids, vitamins and their cofactors, and amino acids were changed in C. sinica larvae.
Amino acids are fundamental for synthesizing proteins and phospholipids, energy production, and involved in morphogenetic processes. In this study, lysine degradation, tryptophan metabolism, phenylalanine metabolism, arginine and proline metabolism, valine, leucine, and isoleucine degradation were substantially down-regulated ( Table 1). The reduced metabolism of these amino acids could significantly impair C. sinica growth and development.
The carbohydrate metabolism is essential for cellular energy balance and for the biosynthesis of new cellular building blocks [43]. In this study, 10 carbohydrate pathways were down-regulated, and eight pathways were upregulated in both ESI+ and ESI− modes (Table 1). Among them, succinic acid was a differentially enriched metabolite in tricarboxylic acid (TCA) cycles. The TCA cycle, known as the citric acid cycle, has an important function that involves the intermediate compounds for the synthesis of amino acids and fatty acids and the formation of large quantities of adenosine triphosphate (ATP) that provides energy for various biological processes [24]. The downregulation of succinic acids of carbohydrate metabolites could cause a shortage of intermediate compounds and energy in azadirachtin-ingested C. sinica larvae, impairing their growth and development. Additionally, the amino sugar and nucleotide sugar metabolism pathway was enriched in azadirachtin-ingested C. sinica larvae. D-glycerate 3-phosphate was also among the differentially enriched metabolites in glycolysis/gluconeogenesis pathway. Glycolysis and gluconeogenesis are metabolic processes responsible for glucose degradation or glucose synthesis, respectively [44].
Ingestion of azadirachtin-treated P. xylostella larvae also affected the metabolism of vitamins and their cofactor in C. sinica larvae. Biotin is an essential substance for insects and affects the development of advanced larvae and pupae. The biotin metabolism pathway was significantly enriched in C. sinica larvae. The downregulation of this metabolite could adversely affect normal metabolism and C. sinica development.
Glycolysis is the main metabolic pathway of carbohydrates, such as galactose and fructose [44]. As an intermediate in both glycolysis and gluconeogenesis, the change in the relative content of D-glycerate 3-phosphate could also affect the generation of intermediate compounds and energy to maintain normal biological processes. Furthermore, the pathways of pentose and glucuronate interconversions, C-5 branched dibasic acid metabolism, ascorbate, and aldarate metabolism were enriched in azadirachtin-ingested C. sinica. Such a series of changes in carbohydrate metabolisms would affect the energy supply of C. sinica larvae, thus their development. Lipids are of vital importance to insects as energy sources and substrates for embryogenesis and development, pupation, metamorphosis, and other activities [45]. They are important components of insect cell membrane and also precursors of many insect pheromones [45]. Azadirachtin could influence the quantity and relative composition of fatty acids [30]. In the present study, 11 lipids or lipid-like metabolites were found to be differentially abundant, of which seven were upregulated, and four were down-regulated. The downregulated included palmitoleic acid, linolenic acid, 9-hydeoxyoctadecadienoic acid (9-OxoODE), and arachidic acid. Linolenic acid plays an important role in insect reproduction as it is a key constitute of oocyte dry mass and the major energy source for embryo development. The linolenic acid has been documented to be required for developing Heliothines subflexa [46]. The reduction in metabolism of linolenic acid in azadirachtin-ingested C. sinica larvae may potentially affect the pupation and eclosion, which required further confirmation.

Key Metabolic Pathways Affected by Feeding on Azadirachtin-Treated P. xylostella Larvae
The aforementioned analyses provide an overall spectrum of metabolite changes in C. sinica after ingestion of azadirachtin-treated P. xylostella larvae. The next question would be which metabolic pathways might be specifically implicated in the reduced percentages in pupation and eclosion C. sinica. Through the KEGG pathway analysis of the differentially abundant metabolites, the disturbed metabolic pathways caused by the consumption of azadirachtin-treated P. xylostella larvae were analyzed. A working model was constructed using the reference map deposited in the KEGG database (Figure 7). It was noticed that 20 differential metabolites were related to amino acid metabolic pathways. Among them, 15 metabolites were down-regulated. These results could imply that the ingestion of azadirachtin-treated P. xylostella larvae might result in the impairment in hydrolyzing proteins in C. sinica resulting in an insufficient supply of amino acids or directly affect amino acid metabolisms. Amino acids, particularly the essential ones are fundamental for insect growth and development [47]. A distinct biochemical characteristic of insects is their higher levels of free amino acids in the hemolymph [48], and the likely utilization of the free amino acids as silk protein synthesis to produce cocoon [49]. Furthermore, the eclosion is controlled by three peptide hormones: eclosion hormone, ecdysis-triggering hormone, and crustacean cardioactive peptide [50,51]. The reduction of free amino acids could hamper both pupation and eclosion of C. sinica.

Precautions When Azadirachtin and C. sinica Are Used in IPM
The changes in a wide range of metabolites in C. sinica larvae suggest that the action mode of azadirachtin is multifaceted with multiple biological targets. Thus, precautions should be taken when azadirachtin and C. sinica are to be used in IPM. Although this study was mainly focused on metabolites without further analysis of the biological influence of these molecular alterations on natural enemies, our data did document that C. sinica larvae after ingestion of larvae of P. xylostella treated with azadirachtin at 2.00 mg/L significantly reduced pupation and eclosion. Such adverse effects were related to the substantial changes in metabolomic profile of C. sinica. As mentioned previously, the premise for the use of pesticides, even botanical ones, in IPM is ascertaining their compatibility with beneficial predators. Our study showed azadirachtin and C. sinica are not compatible, which primarily agrees with the results of C. claveri in responses to azadirachtin [18][19][20][21]. Commercially, azadirachtin has been applied at much higher concentrations ranging from 5 to 100 mg/L [56]. Thus, detrimental effects to beneficial insects could be even severer. However, this does not exclude the use of azadirachtin and C. sinica at different time periods. Nevertheless, our study has raised the compatibility question between the two control tactics. Further studies are needed to identify specific mechanisms underlying the reduced pupation and eclosion of C. sinica larvae, effects on mating and oviposition of adults, and the relationship between azadirachtin concentrations and metabolite changes As shown in Figure 7, a majority of essential amino acid metabolisms were downregulated. A limited supply of tryptophan resulted in the decrease in the contents of indole, L-kynurenine, and 5-hydroxyindoleacetic acid. The downregulation of these compounds could limit the biosynthesis of acetyl-CoA, affecting the TCA cycle. It is worth mentioning that metabolites of the kynurenic pathway generally reach to peak concentrations in insects during pupation [52,53]. The reduced availability of kynurenine could impair pupation in C. sinica. Limited availability of lysine caused reduced biosynthesis of cadaverine and L-pipecolic acid. The downregulation of cadaverine also affected acetyl-CoA biosynthesis, subsequently affecting the TCA cycle. In the arginine metabolism pathway, low arginine in cells resulted in reduced synthesis of N-succinyl-L-glutamate, L-citrulline, and agmatine. Arginine can be converted to proline. Proline was a major substrate used in insect flight metabolism, which is known as the fuel of insect flight [54]. Choline was downregulated, which could indirectly affect the metabolism of serine and threonine. Pyruvate, an indirect derivative of choline, also from glycolysis is regarded as a key metabolite producing valine [55]. Valine also affects acetyl-CoA. L-valine was significantly low in C. sinica. Taken together, at the time of requiring higher levels of free amino acids in the hemolymph, the reduced amino acid metabolisms, particularly the essential ones significantly affected the pupation and eclosion of C. sinica.

Precautions When Azadirachtin and C. sinica Are Used in IPM
The changes in a wide range of metabolites in C. sinica larvae suggest that the action mode of azadirachtin is multifaceted with multiple biological targets. Thus, precautions should be taken when azadirachtin and C. sinica are to be used in IPM. Although this study was mainly focused on metabolites without further analysis of the biological influence of these molecular alterations on natural enemies, our data did document that C. sinica larvae after ingestion of larvae of P. xylostella treated with azadirachtin at 2.00 mg/L significantly reduced pupation and eclosion. Such adverse effects were related to the substantial changes in metabolomic profile of C. sinica. As mentioned previously, the premise for the use of pesticides, even botanical ones, in IPM is ascertaining their compatibility with beneficial predators. Our study showed azadirachtin and C. sinica are not compatible, which primarily agrees with the results of C. claveri in responses to azadirachtin [18][19][20][21]. Commercially, azadirachtin has been applied at much higher concentrations ranging from 5 to 100 mg/L [56]. Thus, detrimental effects to beneficial insects could be even severer. However, this does not exclude the use of azadirachtin and C. sinica at different time periods. Nevertheless, our study has raised the compatibility question between the two control tactics. Further studies are needed to identify specific mechanisms underlying the reduced pupation and eclosion of C. sinica larvae, effects on mating and oviposition of adults, and the relationship between azadirachtin concentrations and metabolite changes as well as threshold concentrations of azadirachtin to C. sinica for potentially better use of the two tactics for pest management.

Chrysopa sinica and Plutella xylostella
Chrysopa sinica and Plutella xylostella were raised in the insect rearing facility of the Key Laboratory of Natural Pesticide and Chemical Biology, Ministry of Education, South China Agricultural University, Guangzhou, China. C. sinica larvae were fed larvae of P. xylostella, and the adults were fed on a diet of P. xylostella larvae, along with 15% honey water, and yeast powder. The larvae of P. xylostella were fed leaves of the host cabbage (Brassica oleracea L.) plants, and the adults were fed 15% honey water. The temperature in the insect rearing facility was 25 ± 1 • C, relative humidity was 60 ± 5%, and a light-dark cycle of 16 h and 8 h.

Chemical Reagents and Instruments
Reagents including methanol, acetonitrile, ammonium acetate, ammonium hydroxide, and formic acid were purchased from CNW Technologies (ANPEL Laboratory Technologies, Inc. Shanghai, China). The internal standards of 2-chloro-L-phenylalanine was purchased from Shanghai Hengbai Biotech Co., Ltd. (Shanghai, China), and azadirachtin (>90%) were provided by Associate Professor Yongqing Tian at the South China Agricultural University.

Experimental Procedures and Samples Collection
Cabbage plants were singly grown in 15-cm pots filled with a peat-based potting substrate. When plants were at a stage of 10 leaves, 20 plants with a uniform growth size were randomly selected. The third instar larvae of P. xylostella were placed on cabbage leaves, 50 larvae per plant. Azadirachtin stock solution (10,000 µg/mL) was made by dissolving it in acetone. The stock solution was diluted with water resulting in a working solution of azadirachtin at 2.00 mg/L, which was sprayed on leaves of 10 plants (5 mL per plant) as azadirachtin treatment (T). The other 10 plants were sprayed with the same concentration of acetone (0.02%) in the same volume per plant as control treatment (CK). After the P. xylostella larvae had eaten the treated cabbage leaves for 12 h (Figure 1a), the third instar larvae of C. sinica were released and fed continuously with P. xylostella larvae from T and CK plants for 5 days, respectively (Figure 1b). The experiment was arranged as a complete randomized design with 10 replications. Five days later, more than fifty larvae of C. sinica were collected from each plant, quickly frozen in liquid nitrogen, and stored at −80 • C. The stored larvae samples were ground into a fine powder in liquid nitrogen and freeze-dried for 24 h until extraction.
To monitor the development of C. sinica, the pupation and eclosion of C. sinica larvae on cabbage plants of the two treatments were observed. After 10 days, when a complete white round cocoon was formed, the pupation was considered complete (Figure 1c). After 20 days, the adults broke out of the cocoon and spread their wings; they were deemed to have completed their emergence (Figure 1d). The experiment was arranged as a randomized complete block design with three replications (three blocks), each block had 50 larvae. The proportion of pupation and eclosion were calculated and analyzed using SPSS software platform (25.0) (IBM Corporation, Somers, NY, USA), and means were separated based on Tukey's HSD (honestly significant different) test at p < 0.01 and p < 0.001 levels.

Metabolites Extraction and Detection
The ground samples of control (CK) and azadirachtin treatment (T), 100 mg each, 10 replicates per treatment, were placed into Eppendorf tubes, respectively. To each tube was added 300 µL methanol and 20 µL 2-chloro-L-phenylalanine, the samples were vortexed for 30 s and sonicated for 5 min in the ice-water bath. The homogenate and sonicate circles were repeated for 3 times, followed by incubation at −20 • C for 1 h and centrifugation at 13,000 rpm and 4 • C for 15 min. The resulting supernatants were transferred to LC-MS vials and stored at −80 • C until analysis.
The quality control (QC) sample was prepared by mixing an equal aliquot of the supernatants from all of the samples to analyze the repeatability of samples under the same processing method. In the process of analysis, one quality control sample was inserted every 6-10 test analysis samples to monitor the repeatability of the analysis process.
The platform for LC-MS analysis consisted of an UHPLC system (1290, Agilent Technologies) with a UPLC HSS T3 column coupled to Q Exactive Orbitrap mass spectrometer (QEO MS) (Thermo Fisher Scientific). The supernatant (200 µL) was taken into the sample bottle (2 mL), respectively.
Mobile phase conditions were set as follows: the mobile phase A was 0.1% formic acid in water for positive, and 5 mmol/L ammonium acetate in water for negative (adjusted the PH value to 9.0 with ammonia), and the mobile phase B was acetonitrile. The elution gradient of the mobile phase was shown in supplementary file, Table S1.
Mass spectrometry conditions included the use of QEO MS to collect MS and MS/MS data with the electrospray (ESI) source conditions set as follows: spray voltage as 3800 V for positive or −3100 V for negative, capillary temperature 320 • C, sheath gas flow rate as 45 Arb, aux gas flow rate as 15Arb, full MS resolution as 70,000, MS/MS resolution as 17,500, strength of collision energy as 3, collision energy as 20/40/60 eV, scanning scope as 70-1000 m/z, scan rate as 7 Hz.

Data Preprocessing and Multivariate Statistical Analysis
The original LC-MS data files were converted to the mzML format by using ProteoWizard and processed by R package XCMS (version 3.2), including retention time alignment, peak detection, peak matching, and peak integration. The data were then filtered by the following criterion: sample numbers containing a metabolite were less than 50% of all sample numbers in a group. OSI-SMMS (version 1.0, Dalian Chem Data Solution Information Technology Co. Ltd, Dalian, China) was used for peak annotation after data processing with in-house MS/MS database.
To initially visualize the differences between different groups of samples, the principal component analysis (PCA) was applied. PCA analysis is an unsupervised multidimensional statistical analysis method that describes the characteristics of the original data set by compressing the original data into countless principal components, which can reflect the overall metabolic difference between each group of samples and the size of variation between the group samples. Partial least squares discriminant analysis (PLS-DA), as a supervised multivariate statistical analysis method, was used to distinguish the metabolomics profile of two groups by screening variables correlated to class memberships in which class memberships were coded in matrix form into Y [57]. Orthogonal projection to latent structures-discriminant analysis (OPLS-DA) is an extension of PLS-DA which incorporates an orthogonal signal correction (OSC) filter into a PLS model. The model quality was assessed based on cross-validation and permutation test [58]. The variable importance in projection (VIP) score of OPLS-DA model and t-test as a univariate analysis were applied to rank the metabolites that best distinguished the different groups in this study. Those with VIP ≥ 1 and a p-value of t-test < 0.05 were considered differential metabolites between two groups [59].

KEGG Pathway Analysis
After the metabolites were found, the metabolites were mapped to KEGG metabolic pathways for pathway analysis and enrichment analysis [60]. The main biochemical metabolic pathways and signal transduction pathways in differential metabolites were analyzed in this study. Significantly enriched metabolic pathways or signal transduction pathways in differential metabolites comparing with the whole background were identified through pathway enrichment analysis.

Conclusions
Azadirachtin as a botanical pesticide has been increasingly used for control of insect pests. This study investigated responses of predator C. sinica larvae after ingestion of P. xylostella larvae treated with azadirachtin, mainly at the metabolite levels. No mortality occurred in C. sinica larvae after consuming azadirachtin-treated P. xylostella larvae, but the percentages of pupation and eclosion of C. sinica were significantly reduced. Metabolomic analysis showed that azadirachtin has effects on the metabolism of amino acids, carbohydrates, lipid, cofactors and vitamins of C. sinica larvae. These effects may impair the growth and development of C. sinica, resulting in reduced pupation and eclosion percentages. Our studies for the first time documented substantial metabolite changes in C. sinica larvae after ingestion of azadirachtin-treated P. xylostella larvae and raise a question about the compatibility between azadirachtin and C. sinica in control of insect pests through IPM.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/metabo12020158/s1, Table S1: The elution gradient of the mobile phase in LC-MS analysis.