2.1 GC–MS analysis of the essential oil and its fractions

The C. myrrha EO was chemically characterized through GC–MS analysis, with a total of identified compounds of 96.50% (Table 1). The EO resulted mainly dominated by furanosesquiterpenes (87.37%), followed by sesquiterpene hydrocarbons (8.04%), oxygenated sesquiterpenes (1.06%), and other compounds in minor amounts. Furanoeudesma-1,3-diene was the most abundant furanosesquiterpene (41.40%), followed by curzerene (25.89%), lindestrene (13.09%), and isofuranodiene (2.00%). Among the sesquiterpene hydrocarbons, β-elemene (2.87%) and germacrene B (2.11%) were the major components. The chemical composition of this EO is linear with those reported for EO from Ethiopian plants [15, 16]. Furanoeudesma-1,3-diene was consistently the major compound across all samples, with percentages ranging from 34 to 39% of the total composition, followed by isofuranodiene (20%) [15], and lindestrene (12–14%) [15, 16]. In contrast, the chemical composition differs from that reported by Morteza-Semnani et al. [17], that found curzerene (40%) as the dominant compound. However, the high levels of curzerene in the EO are misleading since, during GC–MS analysis, isofuranodiene is converted into curzerene through the Cope rearrangement [18]. For this reason, the GC–MS quali-/quantitative analysis of isofuranodiene is generally inappropriate.

Table 1 Chemical composition of the essential oil (EO) of Commiphora myrrha

Then, the EO was fractionated into two main fractions, which were subsequently analyzed by GC–MS. Fraction 1, with 94.37% of identified compounds, was primarily characterized by sesquiterpene hydrocarbons, with β-elemene and germacrene B as the main components (29.08 and 15.81%, respectively). δ-Elemene (6.94%), β-selinene (6.33%), α-selinene (6.12%), γ-elemene (5.64%), and germacrene D (4.79%) were also found in minor percentages (Table S2, Supplementary File 1). The GC–MS analysis of fraction 2 confirmed the presence of furanosesquiterpenes contained in the EO, with a total of identified compounds of 99.70%. This fraction was predominantly characterized by furanoeudesma-1,3-diene (44.95%), curzerene (33.90%), and lindestrene (16.67%) (Table S3, Supplementary File 1), with smaller amounts of atractylon and isofuranodiene (1.85 and 1.60%, respectively).

2.2 Purification of furanosesquiterpenes

Fraction 2 (348 mg) was further purified to give 37.7 mg of furanoeudesma-1,3-diene, 13 mg of curzerene, and 8 mg of isofuranodiene. The structures of the furanosesquiterpenes were confirmed by 1H, 13C NMR, MS spectrometry, and IR spectroscopy aligning with published data [25,26,27]. The full characterization of the isolated compounds is reported in Supplementary File 1 along with the NMR spectra for furanoeudesma-1,3-diene (Supplementary File 1, Fig. S2, S3, S4).

2.3 HPLC–DAD quantitative analysis

The main furanosesquiterpenes of C. myrrha EO and fraction 2, namely furanoeudesma-1,3-diene, isofuranodiene, and curzerene, were quantified through HPLC–DAD analysis (Fig. 1) to avoid the thermal degradation of isofuranodiene and to have a reliable quantification of the compounds. In detail, furanoeudesma-1,3-diene was confirmed as the predominant furanosesquiterpene (68.42 g/100 g EO), followed by curzerene (18.03 g/100 g EO), and isofuranodiene (7.40 g/100 g EO). Similar levels were also found in fraction 2, as reported in Table S4, Supplementary File 1.

Fig. 1

HPLC–DAD chromatogram of Commiphora myrrha essential oil (EO)

These results are consistent with those obtained from the GC–MS quali-/quantitative analysis, where furanoeudesma-1,3-diene was identified as the principal furanosesquiterpene. The presence of curzerene, also confirmed by HPLC–DAD analysis, may be attributed to the extraction protocol applied to C. myrrha oleo-gum resin to obtain the EO, which could have led to the formation of high amounts of curzerene through thermal degradation of isofuranodiene. Further details on the quantification of the furanosesquiterpenes are provided in Supplementary File 1 (Table S4).

2.4 Mosquito larvicidal activity

This study aimed to assess the larvicidal potential of C. myrrha EO on four mosquito species of significant health impact and to identify the class of compounds primarily responsible for the activity. Both C. myrrha EO and its furanosesquiterpene fraction proved to be highly effective against all the mosquito species tested (LC50 from 4.42 to 16.81 and from 3.72 to 5.04 μg/mL, respectively) (Table 2).

Table 2 Larvicidal activity of the Commiphora mhyrra essential oil (EO), and its derived products against Anopheles and Aedes mosquito species

While mortality rates increased with higher concentrations of both the furanosesquiterpene fraction and EO, the former proved to be more effective against all the tested species. However, the furanosesquiterpene fraction achieved a higher mortality rate at lower doses if compared with the EO (Fig. S5, Supplementary File 1), and LC10, LC30, LC50, and LC90 values were reached with lower doses (Fig. S6, Supplementary File 1).

Regarding pure compounds, they were tested on Ae. aegypti since this species showed the highest sensitivity to the EO treatment among the four target vectors. Furanoeudesma-1,3-diene resulted to be the most effective (LC50 3.28 µg/mL), followed by isofuranodiene (LC50 5.58 µg/mL), and curzerene (LC50 7.44 µg/mL) (Table 2, Fig. 2a) (GLMM post-hoc with Bonferroni correction—Furanoeudesma-1,3-diene vs. Isofuranodiene: OR = 22.983, SE = 8.058, z = 8.941, p < 0.0001; Furanoeudesma-1,3-diene vs. EO: OR = 14.484, SE = 5.147, z = 7.478, p < 0.0001; Furanosesquiterpenes fraction vs. Furanoeudesma-1,3-diene: OR = 0.129, SE = 0.046, z = -5.655, p < 0.0001; Curzerene vs. Furanoeudesma-1,3-diene: OR = 0.023, SE = 0.008, z = -10.679, p < 0.0001). Concerning the mortality rates of all the tested products, furanoeudesma-1,3-diene reached the highest mortality rate with the lowest doses if compared with the others, and even a slight increase in its concentration led to a significantly faster mortality trend (Fig. 2b).

Fig. 2

a Overall species mortality rate of Aedes aegypti treated with the essential oil (EO), the furanosesquiterpene fraction, furanoeudesma-1,3-diene, isofuranodiene, and curzerene, regardless of the tested concentrations. b Mortality trend of Aedes aegypti by increasing the concentration of the essential oil (EO), the furanosesquiterpene fraction, furanoeudesma-1,3-diene, isofuranodiene, and curzerene. Colored dots indicate the total number of tested individuals clustered around 0 when alive or around 1 when dead.

This is the first study that reports the larvicidal activity of C. myrrha EO and demonstrates the key role of furanosesquiterpenes in its efficacy. Furthermore, all the tested products exhibited strong activity against the treated mosquito larvae, with LC50 values below 20 μg/mL.

The extracts derived from species of the genus Commiphora Jacq. have been previously tested on mosquitoes. Indeed, the extract deriving from C. myrrha resin showed an LC50 of 281.83 μg/mL on Ae. aegypti larvae [28], while C. caudata (Wight & Arn.) Engl. extracts exhibited mild toxicity on Ae. aegypti, An. stephensi, and Culex quinquefasciatus Say (LC50 ranging from 94.76 to 112.85 μg/mL) [29]. Comparable results were also obtained for the EO from C. berryi (Arn.) Engl. (LC50 ranging from 122 to 175 μg/mL) [30]. Conversely, extracts from C. swynnertonii Burtt resin displayed lower LC50 values (ranging from 3.95 to 27.04 μg/mL) on An. gambiae, Cx. quinquefasciatus, and Ae. aegypti [31]. The EO from C. erythraea (Ehrenb.) Engl. showed a higher toxicity against Cx. pipiens L., Cx. restuans Theobald, and Ae. aegypti, with LC50 values ranging from 10.05 to 29.83 μg/mL [32]. It is worthy of notice that the above-mentioned extracts from Commiphora species are chemically different from the EO tested in this work, being the latter more concentrated in bioactive furanosesquiterpenes.

Regarding other botanical products, Pavela [33] demonstrated that, although many studies report the mosquito larvicidal activity of EOs, only some of them (i.e., Cinnamomum microphyllum Ridl., C. mollissimum Hook. F., C. rhyncophyllum (Miq.), Callitris glaucophylla Joy Thomps. & L.P. Johnson, Auxemma glazioviana Taub., Blumea densiflora D.C., and Zanthoxylum oxyphyllum Edgew.) displayed LC50 below 10 μg/mL [33]. Indeed, among them, only the EOs from A. glazioviana (LC50 of 3 μg/mL against Ae. aegypti) [34] and C. glaucophylla (LC50 = 0.7 μg/mL against Ae. aegypti) [35] showed an efficacy like that reported in this work. Overall, our results demonstrate that C. myrrha EO is an effective product against mosquito larvae.

2.5 GC–MS-driven untargeted metabolomic analysis

The untargeted metabolomic analysis allowed the extraction of 565 ions, 172 of which were putatively annotated and belonged to different classes of compounds such as amino acids, organic acids, sugars and sugar alcohols, and fatty acids among others (Supplementary File 2, Table S7). Data were first explored by unsupervised PCA (Fig. 3a) built under the first two PCs, which explained 90.4% of the total variance (82% of the PC1 and 8.4% of the PC2). As reported by the PCA loading plots, the PC1 was mainly dominated by compounds such as spermidine, 6-hydroxy nicotinic acid, xanthine, hypoxanthine and asparagine, whereas the PC2 by 3-hydroxy propionic acid, 4-hydroxybutyric acid, glycolic acid, 3-phosphoglycerate and lactic acid (Supplementary File 2, Table S7).

Fig. 3

a PCA scores plot between the selected PCs. The explained variances are shown in brackets; b OPLS-DA score plot of all metabolite features; c VIP scores derived from the OPLS-DA model; d Correlation network obtained by DSPC algorithm using metabolites the metabolites discriminating treated from untreated larvae, blue nodes represent metabolites, red lines indicate a positive correlation between metabolites. In contrast, the blue lines indicate a negative correlation. N = 5

The good separation highlighted by PCA allowed further data processing through supervised discriminant multivariate analysis such as OPLS-DA, which is specifically designed to identify variables that contribute to the differentiation between these groups, focusing on the variations that are directly related to group discrimination and enhancing the ability of the model to identify and highlight differences between groups [36].

As reported in Fig. 3b, the validated (Supplementary File 2, Table S7) OPLS-DA model confirmed the separation of the control and treated groups, indicating that their differences were statistically reliable (Fig. 3b). Further, the VIP scores analysis allowed us to identify those compounds that could be considered the main drivers of group separation (Fig. 3c). The model identified 106 compounds with a VIP score higher than 1 (Supplementary File 2, Table S7) and the top 30 were reported in Fig. 3c. Among them, metabolites such as xanthine, 2-hydroxyglutaric acid, hypoxanthine, serine, methionine sulfoxide, and asparagine, among others, were characterized by the highest VIP scores (Fig. 3c).

DSPC is a machine-learning algorithm that discovers the connections between many metabolites derived from few samples [37]. Considering that the low number of samples could be a limitation to the significance of our analysis, we used the DSPC to produce an additional validation of our results. Figure 3d shows the graphical representation of the DSPC results applied to the metabolite intensities obtained from control and treated larvae. Nodes correspond to the metabolites, and the edges represent the correlations among them. Red edges indicate positive correlations, and blue edges are negative correlations. Each node representing a metabolite is defined by its degree and betweenness values (Table 3) [37].

Table 3 DSPC network analysis.

The degree indicates the number of connections that a node has with other nodes, while betweenness measures the number of interconnections. Nodes with high degree and betweenness values are likely to be significant hubs [38]. We used a degree threshold greater than 5 and betweenness greater than 35 to identify metabolites characteristic of the lethal phenotype. Network analysis of control and treated larvae reveals that metabolites primarily involved in amino acid, polyamine, organic acid, and sugar metabolism, such as N-acetyl-D-hexosamine, putrescine, serine, inositol, homoserine, and malic acid, among others, play a central role in the mode of action of EOs (Fig. 3c and Table 3). Moreover, the results highlighted that: i) N-acethyl-D-hexosamine was positively correlated with serine, glucose-6-phosphate, cadaverine and malonic acid; ii) putrescine was positively correlated with glucose-6-phosphate, malonic acid and glycine, whereas iii) serine was negatively correlated with glycine, malonic acid and glucose-6-phosphate, and positively correlated with 3-phosphoglycerate, xylulose, homoserine and N-acetyl-hexosamine (Fig. 3d).

After multivariate and DSPC analysis the annotated metabolites were further analyzed through univariate analysis to highlight all the differentially accumulated metabolites between control and treated larvae (Supplementary File 1 Table S7 and Fig. 4). The t-test univariate analysis revealed that the treatment significantly impacted 137 out of 172 metabolites. Specifically, 56 metabolites showed significant accumulation in treated larvae, while 81 were significantly reduced (Supplementary File 2, Table S7). Further analysis was conducted using Fold Change (FC) with a cutoff of 1.5 (Supplementary File 2, Table S7) to focus on the most strongly affected metabolites (Supplementary File 2, Table S7). These results were then integrated using a Volcano plot analysis, applying a t-test p-value of ≤ 0.05 and an FC > 1.5, as depicted in Fig. 4 and detailed in Supplementary File 1, Table S5. The volcano plot allowed the reduction of the number of affected metabolites to 87, of which 42 were significantly reduced and 45 significantly increased (Fig. 4 and Supplementary File 1, Table S5).

Fig. 4

Important features selected by volcano plot with fold change threshold > 1.5 and t-tests threshold ≤ 0.05. The red and violet circles represent features above the threshold. Note both fold changes and P values are Log10 transformed. The further its position away from the (0), the more significant the feature is N = 5

Among the analyzed compounds, there was a significant accumulation of polyamine spermidine, purines xanthine and hypoxanthine, gamma-aminobutyric acid (GABA), etc. In contrast, the levels of certain compounds, such as the amino acids asparagine and serine, and the organic acid glucaric acid, were significantly reduced in response to the treatment (Fig. 4).

Finally, the data were analyzed through pathway analysis to highlight which pathways were significantly and highly impacted by the EO treatment, considering our metabolomic coverage. The results highlighted that 32 pathways were significantly affected in response to EO treatment (Supplementary File 2, Table S7), but only 21 were characterized by an impact higher than 0.2 (Supplementary File 1, Table S6). Among them, the pathway related to the phenylalanine tyrosine and tryptophan biosynthesis was the most affected highlighting an impact of 1. It was followed by several pathways involved in amino acid biosynthesis (i.e., alanine aspartate and glutamate metabolism, glycine serine and threonine metabolism, arginine and proline metabolism, among others) and sugars metabolism (i.e. sucrose metabolism, fructose and mannose metabolism, glycerolipids metabolism, etc.) (Supplementary File 1, Table S6) (Fig. 5).

Fig. 5

Here are the first four most affected pathways identified in the pathway analysis (Table 3), along with the full list of associated metabolites. Metabolites that belong to these pathways but were not identified during the analysis are highlighted in light blue. The colored boxes, ranging from light yellow to red, represent the statistical significance of the identified metabolites (light yellow indicates a P value closer to 0.05, while red indicates a P value greater than 0.001). Box plots show the trend of each metabolite identified as statistically significant by the t-test, with red boxes representing control larvae and green boxes representing treated larvae. Each code within the boxes corresponds to the KEGG code of the metabolites associated with the pathways; a Phenylalanine tyrosine and tryptophan biosynthesis (C00079—L-Phenylalanine, C00082—L-Tyrosine); b Alanine aspartate and glutamate metabolism (C00025—Glutamic acid, C00026—alpha-Ketoglutaric acid, C00042—Succinic acid, C00049—Aspartic acid, C00064—L-Glutamine, C00122—Fumaric acid, C00152—L-Asparagine, C00334—GABA), c Glycine serine and threonine metabolism (C00037—Glycine, C00065—Serine, C00097—L-Cysteine, C00188—Threonine, C00258—Glycerate), d sucrose metabolism (C00031—D-Glucose, C00085—D-Fructose 6-phosphate, C00092—D-Glucose 6-phosphate, C00095—D-Fructose, C00089—Sucrose, C00208—Maltose)

This research presents a comprehensive study examining the metabolic impact of C. myrrha EO treatment on Ae. aegypti larvae through GC–MS-driven untargeted metabolomic analysis. The analysis suggested that the EO triggered a broad spectrum of metabolic responses. In detail, one of the most impacted metabolisms was that of amino acids and significant changes were observed in pathways related to alanine, aspartate, and glutamate, as well as glycine, serine, and threonine metabolism. A similar metabolic response was found in Tribolium castaneum (Herbst) larvae and Drosophila melanogaster (Meigen) [39] adults after the treatment with organophosphorus pesticides and pyrethroids, respectively. The increase in amino acids suggested protein degradation since, as also reported in our work, an increase in free N-acetylamino acids that are only produced post-translationally by protein N-acetylation was observed [39]. Moreover, in our study, it was observed a decrease in tryptophan, a pivotal molecule involved in protein biosynthesis and whose down-accumulation has been observed in response to several insecticide treatments [39, 40]. In addition, Brinzer et al. [39] further demonstrated that tryptophan catabolism, one of the most impacted pathways in our experiment, is crucial in the defense of the organism against insecticides [39]. The degradation of the above-mentioned amino acid usually leads to an accumulation of neurotoxic compounds [39, 41], such as 3-hydroxykynurenine, which has also been observed in our study after the EO treatment. The reduction of fumaric acid and phenylalanine observed in this study is consistent with what was previously observed by Gao et al. [40]. Indeed, they detected a reduction of these metabolites in the moth Spodoptera frugiperda (JE Smith) exposed to the LD90 of the insecticide spinetoram, and they suggested that the downregulation of phenylalanine may reduce fumaric acid levels, slowing the tricarboxylic acid (TCA) cycle and decreasing guanosine triphosphate production under spinetoram-induced stress.

Besides the changes observed in phenylalanine, tyrosine, and tryptophan, which are end products of glycolysis and aromatic amino acids, in treated larvae it was also observed a general reduction in glycine, serine and threonine content whose metabolism is linked to glucose, glycogen, and pyruvate metabolism, all of which are integral to energy metabolism [41,42,43,44]. The potential alteration of energetic metabolism was confirmed by the impact on sucrose metabolism, the pentose phosphate pathway, and the TCA cycle. Their downregulation could be confirmed by the reduction in glycolytic phosphorylated sugars (i.e. glucose-6-phosphate, fructose-6-phosphate), an accumulation of those belonging to the pentose phosphate pathway (i.e. ribulose-5-phosphate), and the activation of purine metabolism.

The latter has also been observed in An. sinensis Wiedemann larvae in response to the treatment with deltamethrin [42] and has been reported to generate oxygen species (ROS), which could result in DNA damage and apoptosis in the affected individuals [42].

Insecticides, particularly neurotoxic compounds like permethrin, have been shown to disrupt normal physiological processes in Anopheles larvae, resulting in oxidative stress characterized by increased ROS levels. This oxidative stress can lead to cellular damage and, ultimately, death of the larvae [63, 64]. The accumulation of ROS is a well-documented response to various environmental stresses, including exposure to toxic compounds, and plays a critical role in the pathophysiology of insecticide-induced mortality [45, 46]. Earlier research demonstrated that the levels of GABA and polyamines can increase in response to oxidative stress also induced by ROS [47,48,49,50]. The relationship between insecticides, ROS accumulation, and the protective roles of osmoprotective compounds like GABA and polyamines, has been investigated by examining the transcriptomic responses of Anopheles larvae to insecticide exposure. RNA sequencing analyses have revealed significant changes in gene expression related to stress response, detoxification, and metabolic processes, indicating that larvae are actively responding to the oxidative stress induced by insecticides [45, 46]. This adaptive response may involve the upregulation of genes associated with the synthesis of osmoprotective compounds, thereby enhancing the ability of larvae to cope with the detrimental effects of ROS.

2.6 Non-target toxicity

Daphnia magna is a microcrustacean employed as a model organism for acute and chronic toxicity evaluation of compounds in aquatic ecosystems. Its large employment is due to its high sensitivity to toxic agents [51, 52]. Table 4 reports the toxicity of C. myrrha EO, which resulted toxic to D. magna with a LC50 of 4.51 µL/L, and a LC90 of 13.09 µL/L.

Table 4 Effects of Commiphora myrrha essential oil (EO) on Daphnia magna

The assessment of non-target toxicity of insecticides is of crucial importance for their real-world application. Herein, although the observed toxicity resulted marked, it is quite lower than that of conventional insecticides. For instance, pyrethroids are encompassed with a high toxicity on many non-target organisms. Indeed, the LC50 on some fish species at 48 h ranged around 1.6–5.13 µg/L for deltamethrin [53] and 0.05–245.7 µg/L for permethrin [54]. Regarding D. magna, the 48 h effective concentrations (EC50) are 0.16 µg/L for deltamethrin [55] and 0.2–0.6 µg/L for permethrin [56]. Furthermore, some natural products commercially employed as insecticides showed toxicity on this non-target species. For instance, the toxic effect of neem oil has been assessed on D. magna and the EC50 (48 h) value was 0.17 mL/L [