Integrated LC–MS and GC–MS-based untargeted metabolomics studies of the effect of azadirachtin on Bactrocera dorsalis larvae

Azadirachtin exhibits excellent bioactivities against several hundred arthropods. However, current knowlege of its biochemical effect on B. dorsalis larvae is not deep enough. In this study, integrated LC-MS and GC-MS-based untargeted metabolomics were used to analyze the changes of endogenous metabolites and the biochemical effects of azadirachtin on B. dorsalis larvae. Azadirachtin has excellent bioactivities against B. dorsalis larvae in this study, leading to a longer developmental duration, lower survival rate, and low pupa weight. The effect of azadirachtin was investigated on a total of 22 and 13 differentially abundant metabolites in the LC–MS and GC–MS-based metabolomics results, are selected respectively. Pathway analysis indicated that 14 differentially enriched metabolic pathways, including seven influential pathways, are worthy of attention. Further integrated key metabolic pathway analysis showed that histidine metabolism, d-glutamine and d-glutamate metabolism, biotin metabolism, ascorbate and aldarate metabolism, pentose and glucuronate interconversions, and alanine, aspartate and glutamate metabolism in B. dorsalis larvae are significantly relevant pathways affected by azadirachtin. Although extrapolating the bioactivity results in this study to the practical project of B. dorsalis pest management in the field has limitations, it was found that azadirachtin has a significant effect on the primary metabolism of B. dorsalis larvae.

Bactrocera dorsalis is a destructive polyphagous and invasive insect pest of tropical and subtropical fruits and vegetables; this oriental fruit fly has been found to attack many types of commercial fruits and a wide variety of agricultural products 1 . Azadirachtin exhibits excellent bioactivities against agricultural, forestry, medical, and veterinary arthropods [2][3][4] . However, studies on the effects of azadirachtin on B. dorsalis are scarce. Azadirachtin is the main active ingredient in neem. It was reported that neem derivatives are ineffective when used as toxic bait for tephritid fruit flies 5 . Several studies reported that neem seed kernel extracts and azadirachtin deters oviposition of B. dorsalis adults 6,7 . Neem leaf dust significantly reduced the longevity and fertility of B. dorsalis adults by blocking ovarian development 8 . Neem extract could effectively reduce the fecundity, fertility, and post-embryonic development of freshly emerged B. dorsalis flies 9 . However, we found no previous studies in the literature on the activity of azadirachtin against the larvae of B. dorsalis.
The biological effects of azadirachtin include impacts on egg-laying behavior, feeding behavior, growth and metamorphosis, reproduction, activity, and histopathology 10 . The mode of action of azadirachtin against lepidopteran insects can be explained, in part, by effects on digestive enzymes, NADPH cytochrome reductase, and cholinesterase 11 . The physiological effects of azadirachtin include direct inhibition of cell division and protein synthesis 12 . The indirect effects of blocking the release of morphogenetic peptide hormones (PTTH and allatostatins) causes disruption of molting hormone from the prothoracic glands and juvenile hormone from the corpora allata 12 . Transcriptomics analysis to investigate growth inhibition in Drosophila larvae after exposure to azadirachtin showed that 28 genes are significantly up or down regulated, with genes involved in starch and sucrose metabolism, defense response, signal transduction, instar larval or pupal development, and Tr, and QC (quality control) samples were within the 95% Hotelling's T-squared ellipse and significantly separated into clusters. That is no outlier was found among these samples. In the GC-MS analysis, the first principle component (PC1) and second principle component (PC2) explained 28.4% and 27% of the total variance of all samples (Fig. 2a). In ESI− mode, the PC1 and PC2 explained 45.6% and 13.8% of the total variance of all samples (Fig. 2b). In ESI+ mode, the PC1 and PC2 explained 35.9% and 23.9% of the total variance (Fig. 2c).
The supervised partial least squares discrimination analysis (PLS-DA) was performed to identify the metabolites responsible for the separation between control and azadirachtin-exposed groups. The CK and Tr groups in these PLS-DA models were inside the 95% Hotelling's T-squared ellipse and showed clear separation (Fig. 3). The 7-folds internal cross validation and 200 times permutation test were further conducted to assess these models' predictive accuracy and statistical significance. In the GC-MS analysis, the parameters of PLS-DA model's predictive accuracy were R 2 X cum = 0.48, R 2 Y cum = 0.985, and Q 2 Y cum = 0.843; with its corresponding statistical significance were R 2 = 0.793 and Q 2 = −0.0593 (Fig. 3a). In ESI− mode, the parameters of PLS-DA model's predictive accuracy were R 2 X cum = 0.576, R 2 Y cum = 0.999, and Q 2 Y cum = 0.96; with its corresponding statistical significance were 0.844 and 0.0944 (Fig. 3b). In ESI+ mode, the parameters of PLS-DA model's predictive accuracy were R 2 X cum = 0.573, R 2 Y cum = 0.992, and Q 2 Y cum = 0.959; with its corresponding statistical significance were 0.859 and 0.107 (Fig. 3c). According to the criteria that if all blue Q 2 values to the left are lower than the original points to the right or if the blue regression line of the Q 2 points intersects the vertical axis at or below zero 27 , these PLS-DA models exhibited a low risk of overfitting. The above results indicated that these PLS-DA models could identify the differentially enriched metabolites between CK and Tr groups. changed metabolites in B. dorsalis larvae between the CK and the Tr groups. The representative GC-MS and LC-MS total ion chromatograms (TICs) of B. dorsalis larvae tissue samples are shown in Fig. 4. The shape and quantity of peaks between the CK and the Tr groups varied. Approximately 7328 and 13,375 metabolite peaks were deconvoluted in ESI− and ESI+ mode of LC-MS. By contrast, 415 metabolite peaks were deconvoluted in GC-MS. These deconvoluted data were further processed through missing value imputations, filtering, and normalization in MetaboAnalyst 4.0 28 . A total of 1979 and 3904 remaining peaks in ESI− and ESI+ modes in LC-MS and 235 remaining peaks in GC-MS were further annotated using references in existing databases. After the by-products in GC-MS and the exogenous compounds in LC-MS were removed, the differentially abundant metabolites were selected according to the VIP values from the PLS-DA model (VIP > 1.2) and the corrected P values from Student's t-test (q value < 0.05). Table 1 shows that 15 of the 22 differentially abundant metabolites were upregulated in the LC-MS analysis. Table 2 illustrates that two of the 13 differentially abundant metabolites were downregulated in the GC-MS analysis. As shown in Tables 1 and 2, ten amino acids and derivatives of the differentially abundant metabolites were the most differentially abundant metabolites, followed by seven carbohydrates, six lipids, six nucleosides, three organic acids, and two vitamins and cofactors.
Metabolic pathway of differentially abundant metabolites. The KEGG pathway analysis of differentially abundant metabolites was performed by MetaboAnalyst 4.0 to identify the disturbed metabolic pathways caused by feeding with the azadirachtin diet. A schematic overview was constructed using the reference map deposited in the KEGG database (Fig. 5). Eighteen differential metabolites of 14 differentially enriched metabolic pathways present in the B. dorsalis larvae are worthy of attention. We summarized these differentially enriched metabolic pathways into amino acids, carbohydrates, nucleosides, and vitamin and cofactor metabolism.
Among these differentially enriched metabolic pathways, seven had pathway impact values exceeding 0.1, which was the threshold for relevance after the pathway enrichment and topology analysis. On the basis of the negative log(P) and impact values, we characterized histidine metabolism, d-glutamine and d-glutamate metabolism, biotin metabolism, ascorbate and aldarate metabolism, pentose and glucuronate interconversions, and alanine, aspartate, and glutamate metabolism in B. dorsalis larvae as the significantly relevant pathways affected by azadirachtin; their impact values were 0.5, 0.33, 0.25, 0.25, 0.22, and 0.17, respectively (Fig. 6).

Discussion
In this study, azadirachtin mixed in an artificial diet was found to significantly prolong the developmental duration of larvae and decrease the larval survival rate and pupal weight. Although this situation is quite different from practical application in the field because B. dorsalis larvae are inside the fruit, and there is lack of evidence that azadirachtin penetrate the skin and flesh, the metabolomics analysis provides new insights into the biochemical response of B. dorsalis larvae to azadirachtin.
Azadirachtin could significantly reduce the quantity and relative composition of fatty acids 20 . However, only six lipids or lipid-like molecules were found to be differentially abundant in this study, and none of the lipid www.nature.com/scientificreports www.nature.com/scientificreports/ metabolism pathways were enriched. Hence, the targeted metabolomics for fatty acids or lipidomics should be considered to explore a more comprehensive effect of azadirachtin on the lipids of B. dorsalis larvae.
Azadirachtin was found to affect the carbohydrates metabolism of B. dorsalis larvae. Succinic and malic acids were the two differentially enriched metabolites in TCA cycles in this study. The TCA cycle, known as the citric acid cycle, has two important functions. The first involves the intermediate compounds for the synthesis of amino and fatty acids. The other involves the formation of large quantities of ATP, which provides energy for various synthetic processes 29 . The downregulation of such metabolites in B. dorsalis larvae is expected to indicate that a shortage of intermediate compounds and energy in the B. dorsalis larvae fed with azadirachtin diet. Xylitol can be used in the TCA cycle and after conversion into xylulose and arabitol through pentose and glucuronate interconversions (https://www.genome.jp/dbget-bin/www_bget?pathway:dme00040). So changes in the relative content of these sugars is expected to impact the generation of energy and intermediate compounds to maintain normal biological processes.
Azadirachtin was also found to affect the amino acid metabolism of B. dorsalis larvae. Glutamic acid is involved in many biochemical pathways and is regarded as a key metabolite linking carbon and nitrogen metabolism 30 . The relative high content of glutamic acid in treated larvae could be caused by histidine metabolism, d-glutamine and d-glutamate metabolism, and alanine, aspartate, and glutamate metabolism, through which nitrogenous metabolites are further converted into TCA cycles to make up the shortage of energy and intermediate compounds. As histidine regularly plays a key role in the active sites of enzymes 31 , it has to be maintained at a relative high content to meet the needs of enzyme reactions in B. dorsalis larvae fed with azadirachtin diet.
Azadirachtin also affected the vitamins and cofactor metabolism of B. dorsalis larvae, including biotin, also known as vitamin B7 or vitamin H. The relatively high content of biotin in treated larvae in this study could be related to the critical roles this cofactor plays in the intermediate metabolism of gluconeogenesis, fatty acid synthesis, and amino acid catabolism 32 . Myo-inositol is regarded as a vitamin-like essential nutrient 33 . Due to the shortage of intermediate compounds and energy caused by downregulation of succinic and malic acids, myo-inositol was converted into TCA cycles through galactose metabolism, thereby causing its relative low content in this study. Myo-inositol could also be partially responsible for poor growth of B. dorsalis larvae fed with azadirachtin diet, since its deficiency could result in inefficiency in digestion and food utilization and poor growth in shrimps and fish 34 . www.nature.com/scientificreports www.nature.com/scientificreports/ conclusions Although it is inappropriate to extrapolate the bioactivity results in this study to the project of B. dorsalis pest management, the integrated metabolomics analyses have revealed that azadirachtin has a significant effect on the carbohydrates, amino acids, and vitamin and cofactors metabolism of B. dorsalis larvae. This provides new insights into the mode of azadirachtin, the main active ingredient in neem based pesticides.

Chemicals and reagents. Methanol of HPLC grade was obtained from Tianjin Kermel Chemical Reagent
Co., Ltd. (Tianjin, China). Chloroform and acetone were of analytical grade and obtained from the National Pharmaceutical Group (Shanghai, China). The ultrapure water was prepared by ELGA LabWater system.  www.nature.com/scientificreports www.nature.com/scientificreports/ Pyridine, methoxylamine hydrochloride, and N,O-bis(trimethylsilyl)trifluoroacetamide (BSTFA, 98%) with trimethylchlorosilane (TMCS, 1%) were obtained from Aladdin (Shanghai, China). Azadirachtin (>90%) was kindly provided by associate Professor Yong-Qing Tian of the Key Laboratory of Natural Pesticide and Chemical Biology of the Ministry of Education of South China Agricultural University.
Experimental procedures and bioactivities. The population of B. dorsalis was maintained in a laboratory under 25 ± 1 °C, 16:8 h light:dark cycle, and 70-80% RH. The artificial larval diet consisted of corn flour, yeast, sucrose, paper towel, hydrochloric acid, sodium benzoate, banana, and water, whereas the adult diet consisted of water, yeast, and sugar 35,36 .
These experiments contained two groups, namely, control (CK) and 1 μg/g azadirachtin treatment (Tr). The effects of azadirachtin on larval survival, developmental duration, and pupal weight were determined in accordance with a previous report and were briefly summarized below 37 . The stock solution of 10,000 μg/mL azadirachtin was prepared with acetone and stored at 4 °C. For the Tr group, the accurate volume of stock solution was added and fully mixed to ensure that the final azadirachtin concentration in the artificial larval diet was 1 μg/g. The same volume of acetone was added to the CK group. B. dorsalis eggs collected on the same day were counted under a microscope. Three hundred eggs of each replicate were transferred to the diet to investigate larval survival, developmental duration, and pupal weight. Treatments for each group were performed in triplicate.
Sample collection. At the last instar of B. dorsalis larvae (8 days after inoculation), more than 100 larvae were collected and divided into two samples from each replicate of each group. These twelve samples were then snap-frozen in liquid nitrogen and stored at −80 °C. The stored larvae samples were ground into fine powder in liquid nitrogen and freeze-dried for 24 h until extraction.

Metabolite extraction, derivatization, and analysis for GC-MS. The metabolite extraction method
was in reference to a published procedure and described briefly as below 38 . Approximately 50 mg of lyophilized larvae powder from each sample was homogenized for 2 min with 300 µL of precooled solvent (V (chloroform) :V (methanol) : V (water) = 1:2.5:1) for metabolite profiling analysis on GC-MS. The samples were then centrifuged at 12,000 g for 15 min at 4 °C. A total of 250 µL of supernatant was transferred to a new centrifuge tube. The deposit was re-homogenized for 1 min after adding 300 µL of precooled methanol. Subsequently, the samples were centrifuged at 12,000 g for 15 min at 4 °C, and 250 µL of supernatant was incorporated into the first centrifuge tube. The 500 µL combined supernatant was centrifuged at 12,000 g for 15 min at 4 °C. Finally, a total of 450 µL of supernatant was transferred to another new centrifuge tube to make QC and experimental samples.
All samples, including the QC samples, were nitrogen-dried in GC vials at room temperature. The residue was derivatized using a two-step procedure. First, 80 µL of methoxyamine (20 mg/mL in pyridine) was added to the vial and kept at 30 °C for 90 min. Second, 80 µL of BSTFA (1% TMCS) was added to the vial and maintained at 70 °C for 60 min 39 .
The 1 μL of derivatized sample was subjected to Agilent 7890 A/5975 C GC-MS system. It was analyzed in splitless mode with an HP-5MS capillary column (5% phenylmethylsiloxane: 30 m × 250 μm internal diameter, 0.25 μm thickness; Agilent, J & W Scientific, Folsom, CA, USA). The parameters for GC-MS analysis were set in reference to a previous study 40 . First, helium as the carrier gas at a constant flow rate of 1.0 mL/min. Second, the instrument was kept at 50 °C for 1 min, ramped to 100 °C at a rate of 10 °C/min for 1 min, ramped to 200 °C at a rate of 10 °C/min for 1 min, ramped to 280 °C at a rate of 10 °C/min for 1 min, and ramped to 320 °C at a rate of 10 °C/min for 1 min. Third, the injector, transfer line, and ion source were set at 250 °C, 280 °C, and 230 °C, respectively. Finally, mass spectra were acquired with electron ionization mode (70 eV) in full scan mode (m/z 50-650) and solvent delay was set at 5 min. The QC samples were acquired to evaluate stability during analysis. A saturated n-alkane mixture from C7 to C34 was also run at the beginning of the experimental work to determine retention indexes (RIs).

Metabolite extraction and analysis for LC-MS.
Approximately 25 mg of lyophilized larvae powder from each sample was homogenized for 5 min with 800 μL of precooled solvent (V (methanol) :V (water) = 1:1) for metabolite profiling analysis by LC-MS. The samples were then centrifuged at 30,000 g for 20 min at 4 °C. Approximately 650 μL of supernatant was transferred to a new 1.5 mL polypropylene tube and then centrifuged at 25,000 g for 20 min at 4 °C. Next, 550 μL of supernatant was transferred into another new 1.5 mL polypropylene tube to make a QC sample, which was prepared by pooling the same volume of supernatant from each of the samples.
The parameters for LC-MS analysis were set with reference to literature, with some modifications 41 . The samples were subjected to a 2777 C UPLC (Waters, UK) coupled with an Xevo G2 XS QTOF high-resolution tandem mass spectrometer (Waters, UK). The injection volume was 10 μL. Seperation was performed on a HSS T3 C18 column (100 mm × 2.1 mm × 1.8 μm, Waters). The column oven was maintained at 50 °C. The mobile phase consisted of solvents A (water + 0.1% formic acid) and B (acetonitrile + 0.1% formic acid), with a flow rate of 0.4 mL/min. The gradient elution program was set as follows: 0-2 min, 100% A; 2-11 min, 0% to 100% B; 11-13 min, 100% B; and 13-15 min, 0% to 100% A. The Q-TOF mass spectrometer was operated in both positive and negative ion modes. For positive ion mode (ESI+), the detection parameters were set as capillary voltage: 3.0 kV; sampling cone voltage: 40 V; ESI source temperature: 120 °C; desolvation temperature: 450 °C; desolvation gas: 800 L/h; cone gas: 50 L/h; source offset: 80; TOF acquisition mode: sensitivity; acquisition method: continuum MS E mode; TOF mass range: 50-1200 Da; scan time: 0.2 s; and collision energy function 2: trap CE ramp 20-40 eV. For negative ion mode (ESI−), the capillary voltage and desolvation temperature were set at 2.0 kV and 350 °C, respectively. The other parameters were set the same as those in the positive ion mode. During the acquisition, the leucine enkephalin signal was acquired every 3 s to calibrate the mass accuracy. Furthermore, QC samples were acquired to evaluate the stability of the LC-MS analysis during the whole acquisition. www.nature.com/scientificreports www.nature.com/scientificreports/ Data pre-processing and multivariate pattern recognition. The original GC-MS data were automatically analyzed using the automatic mass spectral deconvolution and identification system (AMDIS) and identified by comparing to the database of Fiehn, Golm, and NIST14 42,43 . The RI and mass similarity match were considered in metabolite identification. If the total match factor was greater than 80, then the metabolite identification was reliable 44,45 . Subsequently, the AMDIS output files were extracted and processed using MET-IDEA to obtain the final data set, which included sample information, metabolites, and their intensities 46 . Peaks caused by column bleeding and the BSTFA derivatization procedure were removed.
The original LC-MS data were converted into.abf files by using ABF converter software (https://www.reifycs. com/AbfConverter/). MS DIAL software (http://prime.psc.riken.jp/Metabolomics_Software/MS-DIAL/index. html) with the LC-MS/MS spectral database from MassBank of North America (http://mona.fiehnlab.ucdavis. edu/) was used for peak exaction, data baseline filtering, baseline calibration, peak alignment, deconvolution analysis, peak identification, and peak integration 47,48 . The metabolites were filtered by removing exogenous compounds, such as drugs or compounds of plant origin, on the basis of their likelihood to be present in biological samples 49 .
Total area normalization was performed to reduce the systematic biases within the experiment 50,51 . Data were log transformed and Pareto scaled for multivariate analysis to remove the offsets and adjust the importance of high and low abundance metabolites to an equal level 52 . Multivariate statistical analysis was performed using SIMCA 14.1 demo (Umea, Sweden). PCA showed the distribution of the original data. Supervised PLS-DA was applied to obtain a high level of group separation and to identify the variables responsible for classification 53 .
The PLS-DA model was validated using sevenfold cross validation. The model quality was assessed based on R 2 and Q 2 scores, and the permutation test was conducted to further validate this model 34 . The PLS-DA model was used with the first principal component of VIP values combined with Student's t-test to determine significantly differentially abundant metabolites between CK and Tr. The q values (adjusted P values), which were raw P values from the t-test adjusted using the Benjamini and Hochberg procedure (BH method), were applied to correct for multiple comparisons 45 . The fold change in each metabolite abundance was calculated by comparing the mean values of the peak areas obtained from Tr and CK.
Pathway analysis. Analysis of metabolic pathways affected by azadirachtin was performed using MetaboAnalyst 4.0 28 . This system is a free web-based tool that uses the high-quality KEGG metabolic pathway database as the backend knowledgebase. The hypergeometric test was used for over representation analysis and the out-degree centrality was used for pathway topology analysis. The significantly affected pathways were selected either by P values from pathway enrichment analysis or by impact values from pathway topology analysis 45,54 . The impact values exceeding 0.1 and the negative log(P) values exceeding 2.0 were set as the thresholds in this study.