A completely random sampling method was implemented for the metabolomic study by selecting different collection sites in Colombian locations [Cundinamarca (CU), Nariño (NA), Putumayo (PU), and Quindío (QU)] situated at different altitudes, ranging from (0–3000 m). The geographical distributions of the collected species are shown in Figure 1. All samples were collected from natural bamboo in the period 2020 to 2022 (more details are provided in Supplementary Table S1).

We collected a total of 111 leaf samples from the upper branches of 12 bamboo species, which were distributed across 40 different altitudes and collected under the same conditions. Harvested bamboo leaves were immediately air-dried at room temperature and then ground into powder for extraction. Based on altitude, we divided the 111 leaf samples into two groups: the low altitude group (∼0–1500 m; n = 40) and high altitude group (∼1500–3000 m; n = 71), as indicated in Table 1. At least three individuals (n = 3) were collected at each collection site, and voucher specimens were deposited in the (HPUJ) Herbarium of the Pontificia Universidad Javeriana.

To extract the plant material (leaves), 300 mg of dried and ground samples were mixed with 10 mL of a solvent containing chloroform, methanol, and water (in a ratio of 5:2.5:2.5 v/v/v). The resulting mixture was then vortexed for 1 min, sonicated for 20 min, and centrifuged at 5000 rpm and 20°C for 10 min. The liquid supernatant was filtered using PTFE syringe filters with a pore size of 0.22 µm (Thermo Scientific, Rockwood, TN), and stored at −80°C until used for metabolomics analysis.

Untargeted metabolomics was conducted using an Agilent Infinity 1260 HPLC system coupled to an Agilent 6545 quadrupole time-of-flight (QTOF) mass spectrometer equipped with electrospray ion source (Waldbronn, Germany). 5 μL of the extracts were injected on a C18 column (Kinetex 100 × 2.1 mm, 2.6 µm) at 30°C. The mobile phase employed in LC analyses was composed of 0.1% formic acid in Milli-Q water (v/v) (Solvent A) and acetonitrile (Solvent B) with a constant flow of 400 μL/min. The gradient was set as follows: 3% for 1 min; 3%–97% B in 15 min; 97% B for 2 min; the column was re-equilibrated for 6 min at the initial conditions. Full-scan MS1 and MS/MS spectra were acquired. Data mass spectra were acquired in negative ionization mode (ESI⁻), in a mass range of m/z 80–1700 Da in data-dependent acquisition (DDA) mode. The QTOF instrument was operated in the 4 GHz (high resolution) mode. The parameters used for data acquisition were set as follows: nitrogen used as nebulizer gas with pressure at 52 psi, a capillary voltage of 3000 V, ion source temperature of 250°C, dry gas flow at 12 L/min, and acquisition rate of one spectrum per second. MS/MS fragmentation was performed using a collision-induced dissociation energy of 20 eV. Throughout the analysis, two reference masses were used for mass correction: m/z 112.9856 [C₂O₂F₃(NH₄)] and m/z 1033.9881(C₁₈H₁₈O₆N₃P₃F₂₄).

To evaluate system performance and reproducibility in sample analysis, multiple QC samples were created by pooling and mixing equal volumes of each extracted sample. To assess the instrument’s robustness, pooled QC samples were injected before the sample analysis until system equilibration was achieved and after every ten randomized sample injections.

The LC-QTOF-MS raw data sets were processed using Agilent MassHunter Workstation Profinder software (B.10.0, Agilent Technologies) to extract molecular features for deconvolution, alignment, and integration. The data were manually inspected to eliminate noise and unrelated ions, and a presence filter was applied. For statistical analysis, features that were present in 100% of the samples for each altitudinal group and had a coefficient of variation (CV) of less than 20% in the QC were selected.

Univariate and multivariate analyses were conducted using MatLab (R2019b, MathWorks, Inc., Natick) and SIMCA 14.0 software (Umetrics, Umeå, Sweden), respectively. The multivariate analysis generated PCA and OPLS-DA models, which were validated using cross-validation less than 0.05 and evaluated based on R²X (change in X explained by the model), R²Y (the total of Y explained), and Q² (sum parameter in cross-validation). The significantly differential metabolites were identified by calculating the variable importance in the projection (VIP) greater than 1, with jackknife confident interval (JK) not including zero combined with FC > 2.0 or FC < 0.5. The annotated metabolites and their peak area were organized in a table (.csv), which was uploaded to MetaboAnalyst 5.0 software for statistical, functional, and integrative analysis of metabolomics data (https://www.metaboanalyst.ca/). The software was used for visualization using heatmap clustering (Pang et al., 2022), normalized by Pareto scaling. Additionally, univariate analysis was performed to determine the p-value features.

Differential metabolite annotation was conducted by considering the precision of the mass (maximum error of mass 10 ppm), isotopic pattern distribution, and adduct formation, using different public online databases such as METLIN (https://metlin.scripps.edu/), KEGG (https://genome.jp/kegg), HMDB (https://hmdb.ca/), PubChem (https://pubchem.ncbi.nlm.nih.gov/) and ChEBI (https://www.ebi.ac.uk/chebi/) through the CEU Mass Mediator (https://ceumass.eps.uspceu.es/) tool. The identity of the metabolites was confirmed through MS/MS analysis, which included the use of MS-DIAL 4.80 (https://prime.psc.riken.jp/compms/msdial/main.html), MS-FINDER 3.52 (https://prime.psc.riken.jp/compms/msfinder/main.html), CFM-ID 4.0 (https://cfmid.wishartlab.com/) for in silico mass spectral fragmentation, GNPS web platform (https://gnps.ucsd.edu/ProteoSAFe/static/gnps-splash.jsp) and manual interpretation with the Agilent Mass Hunter Qualitative Analysis software (version 10.0). The metabolites were identified according to the metabolomics standards initiative (Schymanski et al., 2014).

A molecular network was created using the online workflow (https://ccms-ucsd.github.io/GNPSDocumentation/) on the GNPS website (https://gnps.ucsd.edu/). The precursor ion mass tolerance was set to 0.02 Da and an MS/MS fragment ion tolerance of 0.02 Da. A network was then created where edges were filtered to have a cosine score above 0.6 and more than four matched peaks. Further, edges between two nodes were kept in the network if and only if each of the nodes appeared in each other’s respective top 10 most similar nodes. Finally, the maximum size of a molecular family was set to 0, and the lowest-scoring edges were removed from molecular families until the molecular family size was below this threshold. The spectra in the network were then searched against GNPS spectral libraries. All matches kept between network spectra and library spectra were required to have a score above 0.6 and at least four matched peaks (Wang et al., 2016). To enhance chemical structural information within the molecular network, information from in silico structure annotations from GNPS Library Search, Network Annotation Propagation, Dereplicator were incorporated into the network using the GNPS MolNetEnhancer workflow (https://ccms-ucsd.github.io/GNPSDocumentation/molnetenhancer/). Chemical class annotations were performed using the ClassyFire chemical ontology (Djoumbou Feunang et al., 2016; Mohimani et al., 2017; Da Silva et al., 2018; Ernst et al., 2019). The attribute table of the generated nodes was visualized in the Cytoscape software to analyze the molecular network. The data used for the analysis of molecular networks were deposited in the MassIVE Public GNPS database (http://massive.ucsd.edu) with the accession number MSV000090298. The workflow used in this study is summarized in Figure 2.

The quality and stability of the instrument were assessed by conducting a principal component analysis (PCA) on the quality control samples. Supplementary Figure S1 shows a clustering of the QC samples, which indicates the analytical platform’s reliability and the data’s validity. To account for altitudinal variation, we performed PCA and orthogonal partial least squares discriminant analysis (OPLS-DA) models on samples classified into two altitude ranges. Based on this classification, we developed an unsupervised PCA model to observe trends between the two altitude groups. As depicted in Figure 3A, the PCA model showed an initial exploration of the data set, evidencing certain trends between the variables of the altitude groups. Additionally, we constructed a supervised OPLS-DA model to differentiate between low and high altitude groups, with appropriate quality parameters, demonstrating complete separation of the groups. (Figure 3B). This allowed us to identify the significant variables that contributed to the separation between the groups, yielding a total of 89 significant variables with values (false discovery rate, FDR) < 0.05, VIP >1, and FC > 2.0 (or <0.5), as identified according to the metabolomics standards initiative (Schymanski et al., 2014). To validate the OPLS-DA model, we conducted permutation analysis, plotting R ² and Q² of 200 permutation tests, as shown in Figure 3C.

Table 2 presents the differential metabolites (89) that belong to various families, such as flavonoids (48%), fatty acids (19%), cinnamic acid derivatives (11%), unknowns (6%), peptides (3%), glycosylated lignans (2%), alkaloids (2%), carboxylic acids (1%), phenols (1%), carbohydrates (1%), steroids (1%), glycosylated stilbenes (1%) and prenolipids (1%). The largest group of differentially expressed metabolites comprises flavonoids which are the most frequently reported metabolites for these species. Many metabolites in this group showed an upward trend at higher altitudes, potentially contributing to the medicinal potential of these species.

To explore the leading chemical classes, we created a molecular network using the GNPS platform, which enabled us to visualize the chemical space of the metabolome in these species. Figure 4A illustrates the molecular network analyzed in the MolNetEnhancer platform, revealing the classes of MF present in the metabolome. We identified seven groups of molecular families, including phenylpropanoids and polyketides, lipids and lipid-like molecules, organic acid derivatives, lignans, neolignans and related compounds, alkaloids, and derivatives, organoheterocyclic compounds, and oxygenated organic compounds. However, due to the complexity of the metabolome in these species and the limitations of the spectral libraries, many of the nodes did not match any spectral reference. Therefore, it is crucial to continue exploring the chemical composition of these species, which is still limited in the literature. Notably, when we created a molecular network of the flavonoid cluster and compared the low vs. high altitude groups of the annotated compounds, we found that vitexin 2″-O-rhamnoside (FC = 1.494; p-value = 1.41E-03), saponarin (FC = 1.458; p-value = 7.22E-04), and vicenin 2 (FC = 1.360; p-value = 2.13E-02) (Supplementary Figures S2A–C) exhibited positive correlations at higher altitudes, confirming the statistical analysis performed previously (Figure 4B).

By conducting hierarchical clustering analysis and generating a heatmap based on the two most relevant groups of annotated metabolites that exhibited significant changes with altitude (flavonoids and cinnamic acid derivatives), we identified two distinct groups that were associated with altitude. We found that each group exhibited different patterns in terms of chemical composition and levels of presence and abundance. Figure 5 depicts the clustering of the two altitudinal ranges evaluated, with the first cluster corresponding to the high altitude samples that exhibited a high accumulation of flavonoids, while the low altitude samples exhibited a high accumulation of cinnamic acid derivatives.

The box plot shows the differential metabolites for the two altitude groups. In Figure 6 we present the fold change of the main metabolites that exhibited a significant change. We found that cinnamic acid derivatives and flavonoids showed a significant change with the low and high altitude groups, respectively. Our study also highlighted the significance of quinic acid (QA) is a major differential metabolite, which acts as a crucial precursor in the biosynthetic pathway of phenylpropanoids and is essential in the production of a diverse array of phenolic compounds.

To complement the study on the chemical composition of the global metabolome of bamboo species, a molecular network was constructed using the GNPS platform to compare other types of flavonoids that were not statistically significant among the twelve species studied. Four flavonoid C-glycosides (isorhamnetin 7-rhamnoside, isovitexin, isoschaftoside, and rhoifolin) and one flavonoid O-glycoside (cyanidin 3-O-sophoroside) were annotated (Supplementary Figures S2D–H). It was observed that the Guadua angustifolia species was found to have an abundant profile of flavonoid C-glucosides, specifically of the compound isoschaftoside, which is functionally related to apigenin (Figure 7). In terms of flavonoid variation between species, isovitexin and isorhamnetin 7-rhamnoside were found in most species.