The ethyl alcohol (proof, for molecular biology) used for the extraction procedure was acquired from Sigma-Aldrich (St. Louis, United States). The ethyl acetate (HPLC grade) used for the extraction procedure was obtained by J.T. Baker (J.T. Baker-Avantor, Radnor, United States). Acetonitrile and water, both LC–MS grade, were obtained from Fisher Scientific (Fair Lawn, NJ, United States).

Most plant samples were collected by R.F. Almeida on field expeditions throughout Brazil from 2013 to 2017, or were retrieved from discarded fragmented samples used for DNA extraction in molecular studies (de Almeida et al., 2017, 2018; de Almeida and van den Berg, 2020, 2021). For information regarding all sampled specimens, see Supplementary Table 1. After each collection, the plant materials were dried in a desiccator containing silica at room temperature. Samples were then frozen in liquid nitrogen and grounded in a ball mill. The samples were stored in a freezer at −20°C until the preparation of the extracts.

The authorization for conducting this study was conceded by the National System for Management of Genetic Heritage and Associated Traditional Knowledge (SISGEN), n° A6FDC2E.

The plant sample materials were weighed and extracted with EtOH:H₂O (4:1, v/v) or EtOAc (100%) in a proportion of 20 mg of plant material to 1 mL of extraction solvent. The samples were homogenized in a Qiagen TissueLyzer II (Qiagen, Hilden, Germany) for 5 min at 25 MHz and extracted for additional 30 min at room temperature. The samples were centrifuged (5,000 g) for 15 min, and 300 μL of the supernatants were transferred to a 96-deep-well plate. The solvent was dried in a Labconco CentriVap (United States), and the plates were sealed and stored at −80°C prior to analyses.

The extracts were initially resuspended in 200 μL MeOH:H₂O (4:1) containing sulfachloropyridazine (2 μM) as internal standard [to monitor sample injection during the Ultra High Performance Liquid Chromatography (UHPLC)-tandem Mass Spectrometry (MS/MS) data acquisition], and sonicated for 15 min. The plates were centrifuged for 10 min at 1,300 g, and the supernatants were then transferred to a new 96-well plate for metabolomics analyses.

The analyses were carried out with a Thermo Scientific UltiMate 3000 UHPLC system coupled to a Maxis Impact QTOF mass spectrometer (Bruker Daltonics, Germany), controlled by the Otof Control and Hystar software packages, and equipped with ESI source. The extracts were analyzed using a Kinetex 1.7 μm C18 reversed-phase UHPLC column (50 × 2.1 mm; Phenomenex, Torrance, CA, United States), at 40°C, and an injection volume of 5 μL. The pump system consisted of water (A) and acetonitrile (B), both acidified with formic acid (0.1%, v/v), and the flow rate was set at 0.5 ml/min. The metabolites separation was achieved with 5% solvent B for 1 min, followed by a linear gradient from 5 to 100% in 5 min. The column was washed at 100% solvent B for 2 min, then returned to the initial 5% in 1 min, and the equilibration of the column was achieved for 1 min at 5% solvent B. The mass spectra were acquired in both positive and negative ionization modes, separately, in a mass range of 50–2,000 Da in data-dependent acquisition (DDA) mode. The parameters used for data acquisition were set as follows: nitrogen used as nebulizer gas with pressure at 2 bar, a capillary voltage of 4,200 V, ion source temperature of 200°C, dry gas flow at 9 L/min, and spectra rate acquisition of three spectra/s. The five most intense selected ions per spectrum were fragmented (MS/MS) using ramped collision-induced dissociation energy, ranging from 22 to 50 eV. MS/MS active exclusion was set after five spectra and released after 30 s.

The UHPLC–MS/MS data were deposited in the MassIVE Public GNPS dataset² (MSV000085119) and are publicly available.

The raw data files (.d) were converted to .mzXML format using DataAnalysis software (Bruker) after lock mass correction using hexakis(1H,1H,2H-difluoroethoxy) phosphazene (Synquest Laboratories, Alachua, FL, United States), with m/z 622.029509. The quality of the analyses was evaluated considering the retention time and the m/z of a standard solution containing a mixture of six standards, which was analyzed after the completion of each row in a 96-well plate.

The .mzXML files were processed in MZmine2 (version 2.37.corr17.7_kai_merge2) for positive and negative ionization modes, separately. The parameters used for feature finding were as follows: mass detection (centroid, 1.0E3 and 1.0E1 for MS1 and MS2, respectively); chromatogram builder (minimum time span of 0.01 min, minimum height of 3.0E3, and m/z tolerance of 20 ppm); chromatogram deconvolution (baseline cut-off algorithm: minimum peak height: 1.0E3, peak duration range: 0.01–3 min; and baseline level: 1.0E3) with median m/z center calculation, m/z range for MS2 scan pairing of 0.02 Da and retention time (RT) range for MS2 scan pairing of 0.1 min; isotope peaks grouper (m/z tolerance set at 20 ppm, RT tolerance of 0.1 min, maximum charge of 3, and representative isotope set to most intense), join alignment (m/z tolerance of 20 ppm, weight for m/z and RT of 75 and 25, respectively, and RT tolerance of 0.1 min). A filter was applied in order to keep only the features with MS/MS spectra. This feature list was exported as a feature quantification table (.csv), as a MS2 spectral summary (.mgf), and with the SIRIUS export module (.mgf) for downstream analyses.

To investigate the metabolic profile of the dataset, the processed LC–MS/MS data (.mgf and .csv) were used to create a Feature-Based Molecular Network (FBMN) (Nothias et al., 2020) on the GNPS platform (Wang et al., 2016) with input files containing only the features detected in the hydroethanolic extracts. The data were filtered by removing all MS/MS fragment ions within +/− 17 Da of the precursor ion in order to remove possible residual precursor ions, which can sometimes be observed in MS/MS spectra acquired in QToF equipment. Additionally, MS/MS spectra were window filtered to select only the top six fragment ions in the +/− 50 Da window throughout the spectrum. Both the precursor ion and the MS/MS fragment ion tolerance were set to 0.02 Da. A molecular network was created, in which the edges were filtered to have a cosine score above 0.7 and at least four matched peaks. Similarly, the parameters for the library search (for comparison between the experimental and library spectra) were set to have a score above 0.7 and at least four matched peaks to assist in the metabolites annotation—level three according to the metabolomics standards initiative (MSI; Sumner et al., 2007). The FBMN jobs on GNPS can be found at https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=2c5f11403ac847a298e4d7866a491143 (positive mode) and https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=501c16500476451f978311057266fbdf (negative mode).

The molecular network visualization was performed in Cytoscape (version 3.7.2, Cytoscape Consortium, San Diego, CA, United States; Shannon et al., 2003), in which the nodes correspond to ion features, while the edges between the nodes represent the MS/MS cosine scores calculated. Subnetworks in which the nodes were found in significant abundances in blanks were excluded from the Cytoscape visualization to avoid misinterpretations due to contaminants in the analyses. Sample type information was added to color the nodes as pie charts representing the relative abundance of the features across the samples (colors based on the phylogenetic clades A–J). Node size was scaled relative to the sum of the peak areas obtained in the samples in which the feature was detected. Compounds with the same MS/MS spectra, but with different retention times, were represented as separate nodes, indicating isomers.

A chemical hierarchy analysis (Qemistree; Tripathi et al., 2021) was performed with the metabolites detected in the hydroethanolic extraction protocol. For this, we used the q2-qemistree qiime2 plugin,³ in which the feature quantification table (.csv) and the file obtained from the SIRIUS (Dührkop et al., 2019) export module (.mgf) from MZmine were used as input. Briefly, the Qemistree workflow consists of applying SIRIUS (version 4.8.2) to the .mgf file (containing ion information), generating predicted molecular formulas for each feature. The predicted molecular formulas were reranked using ZODIAC (Ludwig et al., 2020), and the predicted molecular fingerprints were subsequently generated using fragmentation trees via CSI:FingerID (Dührkop et al., 2015). The chemical taxonomy of the predicted metabolite structures was obtained by CANOPUS (superclass, class, and subclass; Dührkop et al., 2021). The Euclidean pairwise distances between the molecular fingerprints were calculated, and the fingerprint vectors were hierarchically clustered to generate a tree representing the structural chemical relationships of this dataset. The tree was then pruned in order to keep only the fingerprints classified up to a superclass level. The trees were visualized interactively in EMPRESS (Cantrell et al., 2021) for data exploration, in which clade information was added as relative abundance stacked barcharts to each feature. The dendrogram obtained can be interactively visualized with the .qzv files found in the Github repository.⁴

The feature table exported from MZmine was used to perform unsupervised analysis using Qiime 2 (version 2020.2; Bolyen et al., 2019) bioinformatics pipeline within a Jupyter notebook. The metabolomic profiles were compared using the Bray–Curtis distance metric for comparing different extraction protocols, and using Canberra metric to investigate the four subsets individually (two extraction solvents, and two ionization modes). The Principal Coordinates Analysis (PCoA) plots showing the top three principal coordinates were visualized using EMPeror (Vázquez-Baeza et al., 2013). Permutational multivariate ANOVA (PERMANOVA; Anderson, 2001) was also performed in Qiime 2 on metabolite distance matrices to test for clustering significance (with 999 permutations), and the F statistic was reported as a measure of effect size.

The in silico classes retrieved from the Qemistree workflow (Tripathi et al., 2021) and CANOPUS (Dührkop et al., 2021) were tested for differential enrichment in the most sampled clades in study (clades A, F–J). Initially, the relative abundances were summed for the class group, the data were normalized by arc-sine square root transformation, and the effect of the clades was tested using a simple ANOVA. The adjusted values of p from ANOVA were obtained from the Benjamini-Hochberg method. The pairwise differences between the clades were tested through a post-hoc Tukey test, and the magnitude of differential enrichment was calculated through the log₂ fold change in mean relative abundance between clades. The subclasses statistically enriched for specific clades were then selected to build a heatmap.

In total, 39 genera (out of 74) and 139 species (out of 1,300) of Malpighiaceae were sampled, representing all of its 10 major clades recognized by recent molecular phylogenetic studies (Supplementary Table 1). Since our sampling focused on the diversity of genera and clades within the Malpighiaceae family, genera not sampled on the chemical analyses were not included in the molecular phylogenetic analysis. Sequences for the genes matK, ndhF, PHYC, and rbcL were retrieved from GenBank,⁵ edited using Geneious (Kearse et al., 2012), and aligned using Muscle (Edgar, 2004), with subsequent adjustments in the preliminary matrices to ensure that the nitrogenous bases were correctly aligned. The complete data matrices are available at TreeBase (https://www.treebase.org/treebase-web/search/studySearch.html, accession number S11008).

All trees were rooted in the Byrsonimoid clade (clade A), which is considered the sister-group of the other Malpighiaceae clades according to Davis and Anderson (2010). Combined analysis of plastid + nuclear regions was performed using Bayesian inference and Maximum Likelihood criteria to reconstruct our phylogenetic hypotheses. Both model-based methods were conducted with a mixed substitution model (GTR + G + I) and unlinked parameters, using MrBayes 3.1.2 (Ronquist and Huelsenbeck, 2003) and raxmlGUI2 (Edler et al., 2021). For the Bayesian inference, the Markov Chain Monte Carlo (MCMC) was run using two simultaneous independent runs with four chains each (one cold and three heated), saving one tree every 1,000 generations, for a total of 10 million generations. We excluded 20% of retained trees as “burn in”, and checked for a stationary phase of likelihood, checking for ESS values higher than 200 for all parameters on Tracer 1.6 (Rambaut et al., 2014). The posterior probabilities (PP) of clades were based on the majority rule consensus, using the stored trees, and calculated with MrBayes 3.1.2 (Ronquist and Huelsenbeck, 2003).

Chemical profiles were selected based on the clustering trends of the PCoA analysis for EtOH80 and EtOAc100% extracts, in both positive and negative ionization modes (four subsets). The subsets in which samples clustered into two groups (A and B) represented different chemical profiles. Character coding followed the recommendations of (Sereno, 2007) for morphological phylogenies. Primary homology hypotheses (de Pinna, 1991) were proposed for different chemical profiles, with a total of two subgroups of classes scored (i.e., presence and absence). All characters were optimized on the concatenated tree from the Bayesian inference using the maximum likelihood function with Mesquite 2.73 (Maddison and Maddison, 2006) using the mk1 algorithm.

After annotating the main classes of compounds present in the hydroethanolic extracts, we transformed the obtained quantitative matrices for the positive and negative ionization modes into qualitative matrices considering only the presence and absence of all annotated classes. Since we do not have all the analyzed species of this study sampled in the latest generic molecular phylogeny of Malpighiaceae, we summarized our results at the generic level using an arithmetic mean equation. After summarizing both matrices (retrieved from the positive and negative ionization modes), we compared all classes of secondary metabolites identified in each of them, merging their binary coding into a single column. By doing so, we obtained a total of 78 classes of secondary metabolites that were optimized into a summarized version of the phylogenetic tree presented by us in the previous step, including only the genera sampled in this study. Ancestral character reconstructions were performed using Maximum Likelihood approaches into the software Mesquite 2.73 (Maddison and Maddison, 2006) and visualized in the software Winclada (Nixon, 1999) using the fast optimization method, which favors homoplasies in the analyses.

The extracts (EtOH:H₂O, 4:1 v/v; EtOAc) obtained from the 197 Malpighiaceae samples were analyzed by LC–MS/MS in both positive and negative ionization modes, and the LC–MS data were processed in MZmine2. The feature finding step resulted in a total of 24,440 and 20,150 features detected in positive and negative ionization modes, respectively. Usually, most of the dereplication tools based on LC–MS/MS consider the fragmentation of [M + H]⁺ or [M − H]⁻ adducts. Generally, molecules with pH greater than 7 (basic compounds) could be easily ionized in the positive mode making adducts with proton(s). The formation of deprotonated molecules is usually limited to compounds able to form acidic protons in the negative ionization mode. Although it is of high importance for the characterization of secondary metabolites, the signal intensity for data acquired in the negative mode is usually lower compared to the positive ionization mode. In addition, the use of positively charged ion fragmentation is more relevant due to the larger spectral library availability (Steckel and Schlosser, 2019; Wolfender et al., 2019).

The Venn diagrams obtained in positive and negative ionization modes for the different extraction solvents (Figure 3A) showed that, in the positive ionization mode, 74% of the features were shared between the two extraction protocols, while in the negative ionization mode, this number dropped to 59%, with the number of metabolites exclusively detected in the hydroethanolic extraction rising from 20 to 36%. Only about 5% of the features were detected exclusively with ethyl acetate as extraction solvent, regardless of the ionization mode used. These results are in accordance with previous reports since even though there is a common core metabolome, solvent-specific metabolites are likely to be observed (Crüsemann et al., 2017).

Usually, distinct classes of metabolites are extracted with solvents of different polarities (usually hydroalcoholic mixtures, methanol, ethyl acetate, or methylene chloride), which will allow the enrichment of specific classes of metabolites (such as flavonoids, coumarins, glucosides, alkaloids, diterpenes, and saponins, among others), depending on the solvent polarity (Wolfender et al., 2019; Pilon et al., 2020). In this way, different extraction solvents can be used to obtain a broader chemical diversity. Our results show that the diversity of solvent-specific metabolites will also vary depending on the ionization mode employed.

To observe the chemical space provided by the metabolomic profiles obtained by the different extraction protocols, PCoA of the positive and negative subsets were created (Figure 3B). In summary, differently than Principal Component Analysis (PCA) which measures correlations among the samples, PCoA analysis is used to calculate distances among them, and the way these distances are calculated can result in different clustering trends in the plots. When the Euclidean distance is used in PCoA analysis, the result will be the same as if PCA was employed (Mohammadi and Prasanna, 2003; Bauermeister et al., 2022). To evaluate the impact of different extraction protocols, PCoA plots obtained by Bray–Curtis distance metric showed that the two extraction solvents resulted in very different metabolomic profiles on both ionization modes (for positive ionization mode: PERMANOVA F = 40.39, p = 0.001; for negative ionization mode: PERMANOVA F = 38.37, p = 0.001). Even though most of the features detected are shared between the two extraction protocols, the compounds’ relative abundance may significantly vary resulting in different metabolomic profiles.

In order to determine whether there would be clustering trends for each of the extracts obtained and in both ionization modes, PCoAs for these four subsets were also obtained. For the hydroethanolic extract, two main groups were formed in both ionization modes (Supplementary Figure 1 A; for positive ionization mode: PERMANOVA F = 28.59, p = 0.001; for negative ionization mode: PERMANOVA F = 43.92, p = 0.001). For the ethyl acetate extract, a separation into two groups was obtained in the positive ionization mode (PERMANOVA F = 18.71, p = 0.001; Supplementary Figure 2 A), while no group separation was observed in the negative ionization mode in any of the distance matrices tested (Supplementary Figure 1 B), reinforcing that the extraction solvent and ionization mode influence data acquisition.

The clustering trends observed for the hydroethanolic extracts were optimized using the maximum likelihood criteria into the most recent molecular phylogenetic tree of Malpighiaceae (obtained from Davis and Anderson, 2010), resulting in the cladograms shown in Supplementary Figure 1 B. These trees are graphical representations of how these different chemical profiles obtained from hydroethanolic extracts originated in specific clades (i.e., groups) through time (i.e., geological time). For the metabolites detected in the positive ionization mode, we can infer that the production of metabolites from group B originated in at least three separate geological times (firstly in part of clade A, secondly in clade C, and thirdly in the most recent common ancestor of clade E shared with all remaining clades of Malpighiaceae; i.e., clades E–J). In other words, the samples from group B are, in general, evolutionary more recent than the ones from group A. Some specific genera (such as Amorimia, Alicia, and Banisteriopsis), which are more recent in terms of their origins in past geological times, also cluster together with group A. A possible explanation for this tendency is that these genera are mostly found in dry environments, just like some of the early-diverging genera of Malpighiaceae (Davis et al., 2014; de Almeida et al., 2021).

On the other hand, the metabolomic profile obtained in the negative ionization mode for the hydroethanolic extracts showed a different tendency. The different chemical profiles obtained from group B in the negative ionization mode were only recovered as homoplastic to clade J (i.e., Bronwenia, Diplopterys, Stigmaphyllon, Sphedamnocarpus, and Peixotoa), besides also independently originated in some genera from clades A, E, F, and H.

We also evaluated the metabolomic profile obtained for the ethyl acetate extracts in both ionization modes (Supplementary Figures 2 A,C). Contrary to what we previously observed, the clustering trends observed in the positive ionization mode did not show any clear correlation with the phylogeny of the Malpighiaceae family (Supplementary Figure 1 B), while the negative ionization mode did not show any clustering trends in the metrics tested. These differences in both extraction protocols could be explained due to the different polarity ranges covered by these two solvents. From these results, we observed that different extraction protocols and ionization modes in MS can significantly impact the results from multivariate analysis and chemosystematics investigations.

Since the chemical profiles obtained for the hydroethanolic extracts showed a promising correlation with the evolution of Malpighiaceae, we made an in-depth investigation of the data obtained in both positive and negative ionization modes. The MS/MS library search performed in GNPS resulted in 1,070 spectral matches for the positive and 1,025 matches for the negative ionization modes, resulting in 4.6 and 5.3% of the detected chemical spaces, respectively. These matches were manually evaluated and compared with the literature, resulting in level 2 or 3 annotations according to the MSI (Sumner et al., 2007). In addition, these compounds were also searched in the Dictionary of Natural Products, LOTUS (Rutz et al., 2021), and Scifinder databases to inspect for previous reports in Malpighiaceae.

In addition to the library searches, the MS/MS data were visualized by molecular networking analysis. The molecular families constructed by such analysis represent the similarity of fragmentation patterns obtained by tandem mass spectrometry (MS/MS) analysis. In other words, structurally similar compounds will present similar chemical stability and functional groups, leading to similar fragmentation patterns (Yang et al., 2013). These molecular families consist of nodes (representing MS/MS spectra) and of edges connecting these nodes (representing the cosine similarity between two nodes, which measure the relatedness in MS/MS spectra; Aron et al., 2020). In this way, depending on the cosine score set in the analysis, the connections between the nodes can be more, or less strict. All the molecular families discussed in this section can be found in the Supplementary Material.

The library matches retrieved from the analyses obtained in the positive ionization mode showed the presence of a high diversity of classes of compounds, including C-glycosylated and O-glycosylated flavonoids, lipids, alkaloids, quinic acid derivatives, amides, triterpenes, iridoids, and lignans (Figure 4). Some of these chemical classes were widely detected in all phylogenetic clades, while others were more specific to particular clades or even genera.

Several molecular families related to flavonoid compounds were obtained, highlighting the difference in the fragmentation pattern in MS/MS spectra among these subclasses with different substituents. For instance, the usual MS/MS spectra obtained for O-glycosylated flavonoids in electrospray ionization consist of the neutral losses of the glycosidic substituents, such as hexosides (162 Da), deoxyhexosides (146 Da), and pentosides (132 Da), eventually reaching the aglycone. On the other hand, C-glycosylated flavonoids present a very distinct fragmentation pathway, with many more fragments observed between the precursor ion and the aglycone. The loss of the glycosidic portion is usually reached by the consecutive losses of water molecules and 120 Da, characteristic of C-glycosylated flavonoids (Mannochio-Russo et al., 2020). Therefore, the distinction of these two molecular families in the molecular networking analysis is expected. In addition, compounds containing two C-glycosidic bounds, even more fragments are expected, and we can observe two molecular families for C-glycosylated flavonoids.

Flavonoids with two C-glycosylated portions were mainly detected in samples from clade I, especially from the genera Amorimia and Mascagnia (Supplementary Figure 3). Library matches from this molecular family included apigenin-di-C-hexoside-pentoside and luteolin-di-C-hexoside, being in accordance with the previous reports for the Amorimia genus (Mannochio-Russo et al., 2020). On the other hand, spectral matches to flavonoids containing only one C-glycosylated portion (Supplementary Figure 4) were less clade-specific, being detected in samples from clades A, G, H, I, and J.

In addition to the C-glycosylated flavonoids, three main molecular families were observed for O-glycosylated ones: (1) flavonoids containing only glycosides, (2) glycosides in addition to phenylpropanoids, and (3) glycosides in addition to galloyl portions. These are also explained by the difference in general MS/MS spectra of each of these groups, in which the flavonoids containing the galloyl portion present a characteristic fragment at m/z 153, and the ones containing phenylpropanoids will show the fragment ions regarding this substituent (m/z 147 for coumaric acid, for instance). Finally, O-glycosylated flavonoids containing only glycosides will only present neutral losses regarding each glycosidic portion.

Several matches to O-glycosylated flavonoids were observed in all Malpighiaceae clades (Supplementary Figure 5), which is expected since this class of compounds is considered ubiquitous to plant species (Buer et al., 2010). This network is primarily composed of quercetin and kaempferol derivatives bound to glycosidic portions. A small cluster in this molecular family corresponds to compounds mostly present in clade J, specifically from Stigmaphyllon species. These compounds were annotated as flavonoids with glucuronide, acetylated, and malonylated hexoside substituents, which have not been reported in the literature for the genus Stigmaphyllon to date.

O-glycosylated flavonoids containing sugars and phenylpropanoid portions as substituents were also observed (Supplementary Figure 6), mainly present in the Ptilochaeta genus (clade E), which is in accordance with previous reports (Mannochio-Russo et al., 2020). In addition, sugars bound to galloyl portions (Supplementary Figure 7) were also detected as flavonoid substituents and widely found in all Malpighiaceae clades.

A network containing mainly non-glycosylated flavonoids (Supplementary Figure 8), widely distributed in all clades, was also observed, in which only the fragments relative to aglycones were detected in the MS/MS spectra. These compounds can represent, indeed, aglycones, or even in-source fragments of the O-glycosylated flavonoids. Additionally, a molecular family with spectral matches to catechin and afzelechin and their derivatives (Supplementary Figure 9) was also formed, containing mainly methoxylated portions and galloyl and sugar substituents.

Quinic acid derivatives were detected in higher amounts in clades H and J for phenylpropanoid substituents, and in clades A, B, and G for gallic acid substituents (Supplementary Figures 10, 11). Usually, the MS/MS patterns will contain an ion relative to the quinic acid moiety (usually with additional neutral loss of water), in addition to the fragment at m/z 153, relative to the galloyl portion, while other characteristic fragments will be observed for the phenylpropanoid derivatives (Mannochio-Russo et al., 2020). Previous reports described the presence of galloylquinic acids in Byrsonima (clade A) species (Fraige et al., 2018; Mannochio-Russo et al., 2020), while quinic acids containing phenylpropanoids substituents have been described for the genus Heteropterys (clade H; Huerta-Reyes et al., 2013; Paula-Freire et al., 2013).

Many lipid-like molecules presented library matches, mainly corresponding to the glycerophospholipids (Supplementary Figure 12), fatty acids, and fatty esters (Supplementary Figure 13) classes, besides jasmonic acid derivatives (Supplementary Figure 14), which were largely distributed in all Malpighiaceae clades. Glycerophospholipids also present key fragments in MS/MS analysis, such as the cleavage of the choline group, which can generate a fragment of the choline group itself (if positively charged) and the fragment of the long-chain portion (Ivanova et al., 2007). Lipids represent an important class of compounds, widely found in plants, with key roles in multiple signaling processes (Mamode Cassim et al., 2019). Along with sugars, such compounds are among the main constituents of the oil glands present in Malpighiaceae leaves and flowers (Possobom et al., 2010). Lipids are mainly produced by oil glands that play an important ecological role in Malpighiaceae, a botanical family which is mainly pollinated by oil collecting bees, and the oldest family characterized by oil-bee pollination (Anderson, 1990; Renner and Schaefer, 2010; Davis et al., 2014).

β-carboline alkaloids and other tryptophan derivatives were detected in samples from many clades (Supplementary Figure 15), including Tetrapterys and Banisteriopsis species, corroborating previous reports (Samoylenko et al., 2010; Queiroz et al., 2014). In general, β-carboline alkaloids present MS/MS spectra with fragments relative to the loss of a hydroxyl group (if present) and to the formation of a four-membered ring (m/z 184 and m/z 160, respectively, harmalol as an example). A small cluster into this molecular family showed the presence of compounds without library matches with higher m/z in samples from clade A and J, mainly in Janusia, Banisteriopsis, and Byrsonima species, which can indicate the presence of glycosylated portions due to the mass differences observed between the nodes.

Networks containing library matches corresponding to isoquinoline, protoberberine, and benzylisoquinoline alkaloids were also observed (Supplementary Figures 16–18, respectively). Isoquinoline alkaloids were mainly observed in samples from clades G, H, and J. Berberine alkaloids usually show a retro Diels-Alder (RDA) reaction and a B-ring cleavage, forming fragments m/z 151 and m/z 178, respectively, in the case of scoulerine (Qing et al., 2020). On the other hand, benzylisoquinoline alkaloids will usually present some key fragments, such as the initial loss of the nitrogen atom as ammonia or as methylamine (in the case of methylated nitrogens), in addition to an “even electron”-type McLafferty rearrangement with a reversed charge distribution, and a fragment relative to the benzyl moiety. Therefore, these key fragmentation pathways will generate ions at m/z 299, m/z 192, and m/z 137, taking reticuline as an example (Schmidt et al., 2005). The molecular family representing protoberberine alkaloids showed nodes corresponding to compounds mostly detected in clades H and J, particularly in Stigmaphyllon and Alicia species. Similarly, benzylisoquinoline alkaloids were also mainly detected in samples from clades H and J, especially in Stigmaphyllon species. These classes of compounds have not been described to date for Malpighiaceae species and represent important traits for chemosystematics studies in Malpighiaceae due to their specificity to the above-mentioned clades, especially for protoberberine and benzylisoquinoline alkaloids.

Spectral matches corresponding to amides and polyamines were also observed in molecular families (as shown in Supplementary Figures 19, 20, respectively). These compounds can generate ions relative to the N-cleavage of the amide bond (Barrère et al., 2014). Amides were mainly detected in the Ptilochaeta genus (clade E) as fatty amides and small dipeptides, while the polyamines were primarily observed in samples from clades E, H, I, and J, being this the first report of these classes of compounds for Malpighiaceae.

Several molecular families related to terpenoids were also observed, which are characterized for presenting several fragments in MS/MS spectra depending on the general skeleton of the molecule, such as RDA, McLafferty rearrangement, water losses, among others (Demarque et al., 2016). One of the largest networks obtained (Supplementary Figure 21) showed library matches to triterpenoids and their precursors. These compounds were also largely distributed among all Malpighiaceae clades, in accordance with previous reports. In fact, triterpenoids are the most described class of compounds from Malpighiaceae species to date, with numerous reports in the literature for Acridocarpus, Byrsonima, and Galphimia genera (Cao et al., 2004; Cardoso Taketa et al., 2004; Aguiar et al., 2005).

It is important to mention at this point that the molecular networking approach highlights the chemical similarity of compounds based on their fragmentation patterns, which may not necessarily reflect their biosynthetic origin. For instance, some diterpenoids are grouped in the same molecular family as the triterpenoids, even though their biosynthesis differs from the plastidial MEP and the cytosolic mevalonate pathway. The connection between them can be observed since these two classes of compounds present similar MS/MS fragmentation patterns. A possible way to separate these compounds into two distinct molecular families would be to significantly increase the threshold for the cosine score similarity (set to 0.7 in this work). However, in the molecular networking analysis, the cosine value is set for the entire dataset, and thus, we used an intermediary value that was adequate for most of the families discussed in this work. Since different chemical classes present different key MS/MS fragments, it is very likely that while separating triterpenoids from diterpenoids, other molecular families can be fragmented into smaller networks, even though they can present significant correlations. In addition, it is important to emphasize that even though we can increase this cosine threshold, it is very likely that the diterpenoids and triterpenoids will not be perfectly separated since these two classes share several key MS/MS fragments.

A separate molecular family relative to ecdysteroids, another class of terpenoids, was also observed (Supplementary Figure 22), and mainly detected in samples from clade H, particularly in Niedenzuella species. These compounds are mainly characterized by the presence of consecutive neutral losses of water molecules, in addition to the hydrocarbon side chain cleavage (Lavrynenko et al., 2013; Mannochio-Russo et al., 2020). The presence of ecdysteroids in Niedenzuella multiglandulosa was recently described in the literature (Mannochio Russo et al., 2020; Mannochio-Russo et al., 2020), being the only reports of this class of compounds in the Malpighiaceae family. Our results showed a wide diversity of ecdysteroids in all six Niedenzuella species sampled, being possible chemical markers for the Niedenzuella genus. This represents about 40% of all the Niedenzuella species reported to date, and future studies with other Niedenzuella species should be performed to corroborate this hypothesis. Moreover, these steroids were also detected in Hiraea, Tetrapterys (including in the recently segregated Glicophyllum species, evolutionary close to Niedenzuella), and Peixotoa genera (clades G, H, and J, respectively).

Iridoids, another relevant class of terpenoids, were also annotated based on library matches, which are shown in Supplementary Figure 23. These compounds present key fragments relative to the neutral losses of the glycosidic portions, in addition to water losses and other possible substituents (Wu et al., 2010). The correspondent networks revealed that these iridoids are mainly present in samples from clade J, especially in Stigmaphyllon species, being in accordance with previous reports (Sainty et al., 1981; Davioud et al., 1985). To date, iridoids have only been reported for Stigmaphyllon genus in Malpighiaceae, and in this study, we observe that these compounds are also present in samples from other clades, in particular, Heteropterys oberdanii species (clade H). Furthermore, a molecular family with a spectral match to a secoiridoid (Supplementary Figure 24) was also observed, with part of the compounds widely distributed among all the clades, and part of them mainly detected in Stigmaphyllon species. Finally, molecular families relative to neolignans and furofuranoid lignans (Supplementary Figure 25) were observed and widely distributed in Malpighiaceae clades, which have not been previously reported in this family to date.

The library matches obtained for the negative ionization mode and the molecular networks formed were also inspected, and these analyses showed mainly the same classes of compounds described for the positive ionization mode. Two molecular families stood out for presenting library matches to classes of compounds that were not observed in the positive ionization mode: proanthocyanidins and lignans (Supplementary Figures 26, 27). Proanthocyanidins dimers and trimers were observed in all Malpighiaceae clades, with characteristic fragmentation patterns (such as RDA), corroborating previous reports (Fraige et al., 2018; Mannochio-Russo et al., 2020). Lignans have already been reported for Tetrapterys mucronata (Queiroz et al., 2014), and here we observe that, in addition to T. mucronata species (clade H), these compounds are also present in other clades.

It is important to emphasize that the compounds discussed above represent only a part of the compounds detected in this study. In fact, about 87% of the molecular families (composed of two nodes or more) obtained in both ionization modes did not show any spectral match, including networks composed mainly or exclusively by nodes representing specific clades/genera. These numbers point to the possibility of undescribed natural products. On the other hand, even considering only the spectral matches, it was possible to obtain important information regarding the classes of compounds produced by these plant species. Another crucial point to be considered is that the sampling in this study comprises a larger number of samples from clades A, H, I, and J, and that these clades presented most of the spectral matches observed. It is possible that if a greater number of samples from other clades are included, more spectral matches can be retrieved to give more insights on other clades as well.

In order to amplify the chemical space from the Malpighiaceae dataset and have additional information about the classes of compounds detected, we used the Qemistree workflow combined with the CANOPUS classification tool for systematic compound class annotation (Dührkop et al., 2021; Tripathi et al., 2021). These in silico classifications consist of level 3 annotations according to the MSI (Sumner et al., 2007). In this way, it was possible to construct a chemical tree based on molecular fingerprints from MS/MS spectra and in silico classification tools. A total of 7,489 and 3,773 fingerprints were generated and classified at a superclass level for the positive and negative ionization modes, respectively. In this way, two chemical hierarchy trees were obtained, as shown in Figure 5A.

From the results obtained, it is evident that the ionization mode employed greatly influences the classifications obtained. At a superclass level, the ones most retrieved in the positive ionization modes were the “organic acids and derivatives,” followed by “benzenoids” and “organoheterocyclic compounds.” In contrast, most of the superclasses retrieved in the negative ionization mode consisted of “lipids and lipid-like molecules,” “organic acids and derivatives,” and “benzenoids” (Supplementary Figure 28 A). At a CANOPUS class level, the “carboxylic acids and derivatives” was the main recovered class in both ionization modes. Ion features classified as “benzene and substituted derivatives,” and “azoles” were also observed several times for the positive ionization mode, while for the negative ionization mode, “organooxygen compounds” and “fatty acyls” were the second and third most abundant classifications.

These results confirm some of the conclusions drawn from the molecular networks, such as the presence of lipids and lipid-like molecules in several Malpighiaceae clades. The molecular networks and library searches in spectral libraries, combined with the in silico approaches based on structural databases allowed us to expand the Malpighiaceae chemical space. Several hypotheses raised from the molecular networks were corroborated with the in silico classifications, giving higher confidence in these results.

The classifications obtained in the different ionization modes are also shown distributed in the chromatographic run (Figure 5B). In addition to being possible to observe differences in classifications between ionization modes, the retention time ranges also vary for specific superclasses. For instance, in the negative ionization mode, organic oxygen compounds, and phenylpropanoids and polyketides elute from near the dead volume to approximately 4 min, while the ion features classified as lipids and lipid-like molecules present higher retention times.

The putative chemical classes retrieved as highly correlated with the most sampled clades (clades A, F, G, H, I, and J; ANOVA p < 0.05) were selected to build a heatmap (Supplementary Figure 29). The normalized distribution pattern of the different classes within the sampled genera showed that specific classes are significantly enriched in determined genera/clades, which corroborates our observations retrieved from the molecular networks. In addition, once again, the results obtained for the positive and negative ionization modes differ.

These results show that both ionization modes result in complementary chemical classifications, which is crucial for comprehensive chemotaxonomic investigations. In fact, the in silico tools used in this study relies on public spectral and structural databases, which are known to be more populated with data on the positive ionization mode (Tripathi et al., 2021). In this way, the results obtained for the negative ionization mode are less extensive compared to the positive ionization mode.

Ancestral character reconstruction analyzes have been increasingly encouraged in natural products studies since the early 2010s (Schmitt and Barker, 2009). It has been used in the chemistry of natural products to investigate chemical evolutionary relationships comprising different organisms, such as plants, fungi, and animals (Lumbsch et al., 2006; Bondoc et al., 2013; Allevato et al., 2019; Coley et al., 2019; Chen et al., 2020; Beaulieu et al., 2021). Phylogenetic methods have proven to be a promising approach to explore the evolution of chemical compounds in a specific genus or family of plants and other living organisms. This analysis consists of optimizing binary character states (presence/absence) into a DNA-based molecular phylogenetic tree by statistically testing using Maximum Likelihood Estimation, and recovering which character states (such as chemical classes) characterizes a given clade or taxonomic group within the molecular phylogenetic tree (i.e., genera or major clades recovered by DNA-based phylogenetic studies in Malpighiaceae). In this way, we can determine the statistical probability of the ancestor of a particular genus or clade to exclusively (i.e., synapomorphy) or non-exclusively (i.e., homoplasy) show the presence or absence of a specific chemical class. These exclusive or non-exclusive chemical classes can be used in future studies to circumscribe the analyzed taxa into a new classification, which additionally considers chemical information.

In the present study, the in silico classifications obtained were used to map the evolutionary history of the classes retrieved in our analyses using the maximum likelihood criteria in the Malpighiaceae DNA-based molecular phylogeny (Figure 6). The classifications obtained in the positive and negative ionization modes were combined to provide a general overview of the total number of chemical classes obtained for the family. Classes retrieved as homoplasies or synapomorphies for the major Malpighiaceae clades (A–J) are described in Table 1, while the ones retrieved for all Malpighiaceae clades and genera are extensively presented in Supplementary Table 2. From the 113 in silico classes retrieved, 35 were present in all genera, such as fatty acyls, flavonoids, glycerolipids, phenols, and purine nucleosides, compounds widely distributed in plants, with a variety of ecological roles in these organisms. Nonetheless, future additional studies sampling the remaining botanical families included in the order Malpighiales are needed to properly evaluate the relevance of these 35 classes of metabolites for the Malpighiaceae family as a whole.

Regarding the 10 major clades of Malpighiaceae, all of them were recovered with at least one homoplastic or synapomorphic class of metabolites supporting them (Table 1; Supplementary Table 2). The 10 major clades in Malpighiaceae were characterized by the presence and absence of 22 and 23 classes of metabolites (i.e., homoplasies), respectively, with the absence of glycerophospholipids being recovered as a synapomorphy of Clade D. It is worth mentioning that even though most classes of secondary metabolites circumscribing all the 10 major clades of Malpighiaceae were recovered as homoplasies, their presence or absence was recorded exclusively for each clade, with reversions (i.e., parallelisms) only being recorded within a few distantly related subclades or genera. The results retrieved from this analysis corroborate some of the conclusions obtained from the molecular networks and from the in silico classifications. For instance, glycerophospholipids, fatty acids, fatty esters, furofuranoid lignans, and prenol lipids classes were widely distributed in our molecular networking analyses for all phylogenetic clades, and these classes were also extensively recovered in the ancestral character reconstructions (i.e., appeared in different geological times). In addition, it also clearly shows which classes were important to circumscribe specific clades/genera.

Our chemotaxonomic approach based on MS/MS analyses of Malpighiaceae plant samples allowed us to obtain a comprehensive overview of the classes of secondary metabolites produced by this taxon. Plant secondary metabolites are known to show patterns of occurrence in certain taxa (Wink, 2003). However, it is important to emphasize that the chemical diversity of a sample is highly influenced by many factors, such as genetic variation and environmental influences (soil nutrients, humidity, herbivory, and ecological interactions, among others; Isah, 2019). Studies have shown that even the same species collected in different biomes (Bueno et al., 2021) or different seasons (Zanatta et al., 2021) can produce different relative amounts of specific metabolites. Therefore, even for a single species, a range of factors can be explored to understand how the metabolites are affected to infer their possible ecological roles. However, it is important to emphasize that the ancestral character reconstructions obtained only focus on qualitative characters (presence and absence) rather than quantitative ones.

Our study provides a starting point for follow-up and systematic evaluation of such factors, and in-depth studies must be conducted to confirm and expand these chemotaxonomic conclusions for both Malpighiaceae and Malpighiales. In fact, there is a big jump from a large (but limited) sampling, as the one we present here, toward a chemotaxonomic investigation of the entire plant family. Nonetheless, as shown by several recent studies mentioned above, using phylogenetic methods with chemistry of natural products data is a promising and revolutionary new line of research that aims to elucidate the evolution of specialized metabolites in living organisms. Therefore, future investigations must be conducted to confirm the hypotheses raised in our study, especially for species that do not present any previous phytochemical study. In addition, our efforts were directed to obtain chemical information at a major clades and generic levels. There might be discrepancies if the ancestral character state reconstructions are obtained at different taxonomic levels, such as intrageneric and species levels.

Finally, the relevance of our evolutionary approach to the study of secondary metabolites can be evidenced if we consider as examples three relevant classes of plant secondary metabolites (Table 1: furanoid lignans, isoflavonoids, and piperidines). The absence of the phenylpropanoids furanoid lignans in clades D, G, and I is regarded as an informative homoplasy in our analyzes since this information can be used to point which clades in the family one must focus on in future studies to search for this class of metabolites. This information can also be used in chemosystematic studies of Malpighiaceae to chemically characterize these clades and help establish a new classification system based on chemical compounds and morphology, following de Almeida et al. (2017). On the other hand, the presence of the phenolic isoflavonoids in clades B and E, and the alkaloid class of piperidines in clades G and I are also informative homoplasies that can be interpreted in the same light as the furanoid lignans. In addition, it is important to emphasize that the in silico classifications obtained rely on structural databases populated with compounds from diverse sources, including plants and microorganisms metabolites. Therefore, some in silico classifications of metabolites more usually found in microorganisms may occur in plant datasets (for instance, the features classified as lactams). A deeper investigation in each taxonomic group should be performed in the future to confirm their presence. In fact, the population of such databases with more compounds derived from plants is necessary to have more accurate information regarding plant species in the future.

The chemical characterization of all Malpighiaceae clades is the first step toward enabling a new research line on the evolution of secondary metabolites in this plant family since this family already has a dated and calibrated molecular phylogeny available in the literature (Wink, 2003; Davis et al., 2014). Merging both analyses would allow us to infer the geological time in which all classes of secondary metabolites have arisen in Malpighiaceae and correlate these dates with past biogeographic events, such as the colonization of different biomes by the most recent common ancestor of all analyzed clades (de Almeida et al., 2018). Additionally, it would also be possible to evaluate which Malpighiaceae lineages experienced a higher diversification throughout the geological time and identify, and which classes of secondary metabolites are correlated with these diversification events (Xi et al., 2012).