Chemical properties of foliar metabolomes represent a key axis of functional trait variation in forests of the tropical Andes

Plants interact with their environment through diverse specialized metabolites that protect them from abiotic stressors like drought or radiation and biotic stressors like herbivores or pathogens. However, few studies have considered the chemical properties of metabolites as a potential axis of functional trait variation along environmental gradients.  Here, we examined how the chemical properties of foliar metabolomes, such as mean aromaticity, hydrophobicity, and polarity, as well as commonly used morphological traits, vary with climate and elevation among 16 forest plots in the tropical Andes of Bolivia. We found that chemical properties were weakly related to morphological traits among tree species, yet both varied significantly with climate and elevation. In particular, abundance-weighted mean hydrophobicity decreased, and polar surface area increased with elevation and in colder and drier climates. Additionally, co-occurring species showed increasing chemical similarity with elevation for the most-aromatic and most-polar metabolites. These results suggest that abiotic stress associated with colder, drier climates and solar radiation acts as a filter for metabolome chemical properties. This contrasts with chemical dissimilarity observed at lower elevations, which is likely driven by pressure from host-specialized enemies in warmer, wetter climates. Our results introduce the possibility that chemical defenses may be constrained by abiotic stressors. Morphological traits and foliar metabolome chemical properties for each species-by-plot are reported in Dataset S1. Community-weighted mean values are reported in Dataset S2. The structural similarities among 20,571 metabolites are reported as a Qemistree dendrogram in .tre phylogeny format as Dataset S3. Masses, molecular formulae, predicted structures, classifications, and chemical properties and sample-level abundances for 20,571 unique metabolites are provided in Dataset S4.

Methods
Forest plot data: The Madidi Project

            Floristic data were collected as part of the Madidi Project (www.mobot.org/madidi), a collaboration of more than two decades between the Herbario Nacional de Bolivia and the Missouri Botanical Garden to document the flora of the Madidi region in the Andes of Bolivia (35). The region features wide variation in plant communities over an extreme elevational gradient, from lowland rainforests located at 200 m above sea level (a.s.l.) to alpine environments above the tree line at 6,000 m a.s.l. (48). The Madidi Project includes 50 1-ha permanent forest plots ranging in elevation from 212 m to 3334 m a.s.l. We selected 16 1-hectare (ha) permanent plots in which leaves were sampled for chemical analysis and which represent broad variation in elevation (662-3324 m a.s.l.), climate, and tree species richness (17-137 species per 1-ha plot). The 16 plots include three seasonally dry, low-elevation forest plots and 13 moist, montane forest plots (28). Abundant genera include: Miconia (Melastomataceae), Sloanea (Elaeocarpaceae), and Ocotea (Lauraceae) in the low-elevation moist plots; Weinmannia (Cunoniaceae), Hedyosmum (Chloranthaceae), and Clethra (Clethraceae) in the high-elevation (>2500 m) plots; and Weinmannia, Hedyosmum, and Clethra in the seasonally dry plots. Tree species richness declines with elevation among the 13 moist forest plots, whereas the three seasonally dry lowland plots display low species richness (28). In each 1-ha plot, all free-standing woody plants with a diameter at breast height of ³ 10 cm were mapped, measured, and identified to a botanically valid species or morphospecies.

Morphological Functional Traits

           Protocols for morphological functional trait data collection are described in detail in the Madidi Project manual (www.mobot.org/madidi). We selected five morphological leaf and stem traits that reflect a species position on a tradeoff axis from conservative traits associated with defense and survival to acquisitive traits associated with fast growth (Box 1). Leaf area and specific leaf area (SLA; area per unit mass) are associated with a resource-acquisitive life history strategy; leaf thickness, bark thickness, and twig specific density are associated with a resource-conservative life history strategy (49). Morphological traits for each species-by-plot are reported in Dataset S1. Community-weighted mean trait values are reported in Dataset S2.

Chemical Analysis

            We collected leaf samples from 473 tree species representing 906 unique species-by-plot. Within each forest plot, we collected leaf samples from 62-90% of the species in the plot (mean = 80% of the species per plot; (28)). Leaves of up to five individual trees per species per plot were collected between 2010 and 2019 and dried with silica gel upon collection in the field. Leaf samples were extracted for untargeted metabolomics analysis following Sedio et al. (31). Briefly, 50 mg of dried leaf tissue was ground to a fine powder and 10 mg weighed for extraction in 1800 mL 90:10 methanol:water pH 5 overnight at 4 °C. Extracts of up to five individuals per species per plot were pooled to create 906 pools representing unique species-by-plot.

All individual extracts and species pools were filtered and analyzed using ultra-high performance liquid chromatography-heated electrospray ionization-tandem mass spectrometry (UHPLC-HESI-MS/MS) using a Thermo Fisher Scientific (Waltham, MA, USA) Vanquish UHPLC with a C18 column and a Thermo QExactive quadrupole-orbitrap MS. Separation of metabolites by UHPLC was followed by HESI ionization in positive mode using full scan MS1 and data-dependent acquisition of MS2. Detailed instrumental methods are described by Sedio et al. (31). Spectra were deposited as a public MassIVE dataset on the Global Natural Products Social (GNPS) Molecular Networking server (doi:10.25345/C52R3P21H).

Raw spectra were centroided and processed for peak detection, peak alignment, and filtering using MZmine 2 (50). Aligned chromatograms were used to create a feature-based molecular network (FBMN; (51)) using GNPS (52). The structural similarities of all metabolites as represented in the resulting network were used to create a dendrogram using the software Qemistree (53), which is reported in Dataset S3. Metabolites were annotated by predicting molecular formulae using Sirius (54), predicting molecular structures using CSI:FingerID (55) and classifying compounds using Canopus (56) according to the organic chemical taxonomy scheme of ClassyFire (57) and according to biosynthetic origins using NPClassifier (58). For a comparison of intra- and inter-specific variation for selected species-rich high- and low-elevation genera, see (28).

To calculate chemical properties of metabolites, we used the highest-confidence molecular structure predicted by CSI:FingerID, represented as a SMILES text string, to query the Chemistry Development Kit (CDK; (59)) using the R package ‘rcdk’ (60). The CDK library includes 51 variables that describe chemical and physical properties, but Walker et al. (27) found that many of these are highly correlated and hence represent five major axes of variation. A correlation matrix of 21 chemical properties for metabolites in our data closely matched that of Walker et al. (27). Hence, like Walker et al. (27), we chose one of each of five major dimensions of variation (Box 1). Molecular formulae, predicted structures, classifications, and chemical properties and sample-level abundances for 20,571 unique metabolites are provided in Dataset S4. Foliar metabolome chemical properties for each species-by-plot are reported in Dataset S1. Community-weighted mean values are reported in Dataset S2.

We calculated the chemical structural-compositional similarity (CSCS) of species, which accounts for the structural similarity of unique metabolites (30). We calculated CSCS with respect to metabolites in the upper and/or lower quartile of nAtomP, ALogP, TopoPSA, and Fsp3, respectively, for the species co-occurring in each of the 16 forest plots.

Climate Data

            We selected four climatic variables to represent variation among the 16 forest plots in temperature, precipitation, and seasonality. Annual mean temperature and annual range in temperature were derived from WorldClim Version 2.1 (61). Annual precipitation and precipitation seasonality, calculated as the ratio of the standard deviation to the mean precipitation of each month, were derived from the Tropical Rainfall Measuring Mission (TRMM), a regional database that provides greater accuracy in precipitation measurements relative to WorldClim in the Bolivian Andes (28). We scaled and centered the four variables and carried out a principal components analysis, of which the first principal component represented 71.2% of the variation and was clearly interpretable as a gradient from cold, dry environments (values < 0) to warm, wet environments (values > 0; (28)). Elevation and position on climate PC1 for each of the 16 forest plots are reported in Dataset S2.

Discipline-Specific Metadata

The DisciplineSpecificMetadata.json file contains parameter values for experimental and instrumental protocols used in liquid chromatography-mass spectrometry (LC-MS) data collection. These methods are also reported in Sedio et al. 2021 "Chemical similarity of co-occurring trees decreases with precipitation and temperature in North American forests". Front. Ecol. Evol. 9.679638. doi: 10.3389/fevo.2021.679638