C. calisaya trees grow as isolated individuals or in clusters of a few trees and due to overharvesting and deforestation, C. calisaya trees can no longer be found at many historical sites (Ruiz, 1792; Andersson, 1998; Honigsbaum, 2001) present. Consequently the trees are increasingly hard to find. We revisited sites in the Yungas region of Bolivia, where most historical records occur (Andersson, 1998; Maldonado et al., 2015, Figure 1), and sampled C. calisaya trees in October 2012 and 2014 (Figure 3; Table 1).

A total of 22 specimens were sampled for this study (Table 1). The sampled trees varied in height from 1 to 7 m with 3 m as average (Table 1) and for each individual, three samples of bark were collected at three different heights of the trunk (0.7, 1.3, and 2.0 m; Figure 2b) for chemical analysis to account for possible variation along the tree (Hodge, 1946; Nair, 2010). Trees below 2 m were only sampled at available heights.

Leaf samples were collected in silica for DNA analyses. Three soil samples were collected around each tree at 10 cm of depth (see Environmental variables section for details). The voucher specimens are deposited in the Bolivian National Herbarium (LPB) and duplicates in the Natural History Museum of Denmark, University of Copenhagen Herbarium (C).

Alkaloids were extracted from all bark samples across the three heights and a test of similarity in content of alkaloids across these three heights was conducted. A few bark samples were covered in lichens, which were brushed off before the bark samples were ground to homogeneity using a custom made grinding device (Hansen et al., 2015), and 50.0 mg of the ground bark was dissolved in 1.0 mL DMSO (dimethyl sulfoxide), and ultrasonicated for 10 min. Following Holmfred et al. (in press), 4.0 mL of 70% methanol containing 0.1% formic acid was added and the mixture was ultrasonicated for 10 min. After centrifugation at 5,000 rpm for 10 min, the clear liquid phase was transferred into a 50.0 mL volumetric flask. The extraction was repeated and the liquid phase transferred into the same volumetric flask and the total contents diluted to 50.0 mL with a solution of 0.1% formic acid in demineralized water (Milli-Q). Recovery from the two extractions was 70–73 and 26–28% respectively (Holmfred et al., in press).

Whereas the official assay of the European Pharmacopoeia for Cinchona bark (Council of Europe, 2016) detects the major alkaloids as two pairs of diastereoisomers, we used an optimized extraction and HPLC method allowing separation and quantification of all four major alkaloids (Figure 2c; Holmfred et al., in press). HPLC analyses were performed with an Agilent 1200 system consisting of an on-line degasser G1379B, a binary pump G1312B, an autosampler G1367C, a column oven G1316B, a DAD detector G1315C, and a FLD detector G1321A. The column used was a KINETEX XB-C18 (150 × 2.1 mm) with 2.6 μm particles. The gradient elution was performed using mobile phase A: 0.2 M ammonium formate buffer with pH 3.0 and water (10:90 v/v) and mobile phase B: 100% methanol at a flow of 0.2 mL min⁻¹. The gradient was 18% B from 0 to 10 min, then changed to 35% B from 10 to 25 min returning to 18% B after 26 min with a total run time of 40 min. The column oven temperature was 50°C and the injection volume 3.0 μL. For UV detection: 250 and 330 nm and fluorescence detection: Excitation 330 nm and emission 420 nm was used. Quinine sulfate and cinchonine standards were obtained from Merck (Darmstadt, Germany). Quinidine and cinchonidine were both obtained from Fluka (Sigma-Aldrich, Denmark). Limit of quantification (LOQ) of the method is 5 μg/g for the four major alkaloids and repeatability of the method, expressed as the relative standard deviation (RSD) of the reference standards, was typically found to be less than 5% (Holmfred et al., in press). Total major alkaloid content was summarized and expressed in mg/g (Table 4).

To assess phylogenetic patterns of plant defenses, we produced a molecular phylogeny including 22 specimens of C. calisaya and three specimens of Cinchona micrantha Ruiz & Pav. as outgroup. Total genomic DNA was extracted from 20 mg dried leaf material using a DNeasy Plant Mini Kit (Qiagen, Denmark) with lysis of the DNA performed by incubating in a standard lysis buffer according to the manufacturer's protocol.

Four chloroplast (trnL-F, matK, rps16, and ndhF) and one nuclear (ITS) DNA regions were sequenced using the primer sets listed in Table 2. Selection of sequence regions was based on previous studies of related taxa showing useful proportions of phylogenetically informative sites within the tribe (Andersson and Antonelli, 2005; Manns and Bremer, 2010). PCR amplification was performed in a 25 μL reaction mixture containing 12.5 μL GMIX (Epicentere an Illumina Company, Madison, WI, USA), 9.5 μL H₂O, 1 μL 20 μM of each primer, 1 μL DNA and 0.125 μL of VWR Taq DNA polymerase (VWR international) using the PCR programs described in Table 2. PCR products were purified using the QIAquick PCR Purification Kit (Qiagen, Denmark) before sequencing on an AB3130 automated sequencer (Applied Biosystems, Foster City, CA, USA). All products were sequenced in both forward and reverse directions using BigDye v 3.1 sequencing reagents and protocols (Applied Biosystems, Foster City, CA, USA). Sequences were assembled using GENEIOUS v R6-7 (www.biomatters.com). All sequence data were submitted to NCBI/GenBank and the accession numbers are listed in Table 3.

Sequences were aligned using MAFFT (Multiple Alignment using Fast Fourier Transform) v 7.017 (Katoh et al., 2002) as implemented in GENEIOUS. Using the same program, the alignment ends were trimmed manually. Data were partitioned into nuclear (ITS) and plastid (trnL-F, matK, rps16, and ndhF) regions and each region was allowed partition-specific parameters. The best suitable evolutionary model was selected according to the Akaike information criterion (AIC) and conducted independently for each partition using the program jModelTest (Guindon and Gascuel, 2003). Phylogenetic relationships were inferred using a standard Bayesian approach as implemented in MrBAYES v 3.2.2 (Huelsenbeck and Ronquist, 2001). Analysis consisted of 3,000,000 generations, two chains, and a sampling frequency of 100 (mcmc NGEN = 3,000,000, SAMPLEFREQ = 100, NCHAINS = 2). Trace files were assessed in TRACER v 1.5 (Rambaut and Drummond, 2007) to ensure stationarity had been obtained and to select an appropriate burn-in phase (in this case the first 10,000 sampled trees).

Phylogenetic signal was assessed for each of the four major alkaloids independently, as well as for cumulative total major alkaloid content using Blomberg's K (Blomberg et al., 2003), a standard index developed for continuous traits. K values of 0 (K = 0) suggest random dispersal of the trait across the tree and therefore no correspondence with phylogeny, while values greater than one (K > 1) suggest a correlation between the phylogeny and the evolution of the trait under Brownian motion. Statistical significance was tested by random shuffling of tips one million times using the 50%-Majority Rule Bayesian consensus phylogeny and the PHYLOSIGNAL function in the PICANTE R package (Kembel et al., 2010; Venables et al., 2014). A 1,000 trees were extracted from within the 95% highest posterior density region of the Bayesian trees, randomly subsampled, and Blomberg's K calculated for each phylogeny using R.

In the present study, “environmental parameters” refer specifically to three types of abiotic macro-environmental variables (Moore et al., 2014): Climate, soil properties, and geographical distance.

Climate variables were taken from WorldClim data (Hijmans et al., 2005). Mean annual temperature and precipitation were downloaded to 30 s resolution based on GPS coordinates of the samples, as measured in the field. Additionally, due to the large altitudinal range of the study species (200–3300 m, Figure 1), altitude and geographical coordinates were recorded by GPS (GARMIN eTrex H) in the field for each sample.

For the soil properties analyses, each of the three subsamples were sieved to under 2 mm before 35 g of the resultant powder was pooled together and sent for analysis by Eurofins (Denmark) for pH, humus content, and content of carbon (total), magnesium, nitrogen (total), phosphorous, and potassium. All methods are standard environmental analyses following the guidelines from the Danish Agrifish Agency under the Ministry of Environment and Food of Denmark (Sørensen and Bülow-Olsen, 1994). Humus was determined with an automated analyzer following the method of Ter Meulen as CO₂ formed by combustion in a muffle furnace (Nabertherm, Germany) in a CO₂ free environment (LECO Tru Mac N, Michigan). Total nitrogen was determined by LECO Tru Mac N as ammoniacal nitrogen NH₄-N by distillation and titration following destruction of organic matter with sulfuric acid. Phosphorus was determined spectrophotometrically using a Foss, FIAstar 5000 Analyzer (Olsen's P), Potassium by flame photometry and Magnesium spectrophotometrically using (AA) ICP-OES, Thermo Scientific, iCAP 6500 Radial, by inductively coupled plasma optical emission spectrometry. Soil pH was determined using an ACCUMET AB15 (Fischer Scientific, UK) and automated with a custom built robot.

The alkaloid variation was analyzed against environmental variables. Due to the large number of explanatory variables compared to dependent variables, principal component analysis was performed on climatic and soil properties (Borcard and Legendre, 2002) to reduce the number of explanatory variables. Additionally a Principal Coordinates Neighborhood Matrix (PCNM) was established from GPS coordinates for each sample. Principal components were used in subsequent multivariate analyses, with components added until over 90% of variation in each ordination was accounted for. Thus, a single component from the climatic ordination, two from the soil and two geographical, were analyzed against the concentration of the major alkaloids with PERMANOVA, using the Adonis function of the vegan package in R and 9,999 permutations (Oksanen et al., 2013). To test for any significant effect of sampling over different years on the alkaloid production of C. calisaya, sampling year was included as an explanatory factor in a PERMANOVA test to account for effects of sampling over different years. Additionally, as much of the variation of individual parameters is lost in performing ordinations used within the PERMANOVA, correlations between individual alkaloids (and total major alkaloid content), and individual climatic parameters using Spearman's rank-order correlation coefficients were used to expand upon the results found within multivariate analysis.

Due to the presence of significant phylogenetic structure, phylogenetic generalized least squares (PGLS; Martins and Hansen, 1997) analyses were performed on the environmental parameters that significantly correlated with C. calisaya alkaloid content, partitioning their effects from those of phylogeny. PGLS was performed using 1,000 trees taken from the 95% posterior density region of the Bayesian trees, randomly subsampled, as with testing for phylogenetic signal with Blomberg's K.

Our sampling sites covered a range in altitude with samples collected from 747 to 2,092 m above sea level (Table 4). Annual precipitation also varied greatly between sampling locations, from 558 to 1,810 mm year⁻¹, and correlated negatively with altitude (r = −0.620, P < 0.001). Despite average annual temperature varying between 10.3 and 17.5°C across sites, temperature did not correlate with altitude or precipitation (r = −0.132, P = 0.269, and r = −0.107, P = 0.309).

Soil properties varied greatly among sites, with humus content varying between 2.5 and 40.0%, C between 60.5 and 97.5%, and N between 0.1 and 1.2%. The soil was highly acidic, as typical of tropical regions, with pH ranging between 4.0 and 6.1. Low concentrations between 0.0 and 1.2 (mg kg⁻¹) were found for P. Both K and Mg varied considerably, K between 3.9 and 26.0 mg kg⁻¹, and Mg between 2.0 and 88.0 mg kg⁻¹ among the sites (Table 4).

The combined aligned DNA matrix comprised 3,700 base pairs (bp) derived from the five DNA regions: ITS (543 bp); trnL-F (839 bp); rps16 (746 bp); matK (764); and ndhF (808 bp). The topology derived from Bayesian analyses is shown in Figure 3 with Bayesian posterior probabilities indicated above the branches. We obtained two major clades of C. calisaya (A and B). Clade A (PP = 1.00) includes trees with high alkaloid content (13–41 mg g⁻¹ cumulative of four major alkaloids) restricted to 1,100–1,350 m altitude between Chulumani and Coripata in Bolivia, except for CMG4214, which was collected near Coroico (Figure 3; Table 4). Clade B (PP = 1.00) primarily includes trees with low alkaloid content, but with a few exceptions. Clade B is a higher altitude clade ranging from about 1,900 to 2,300 m altitude (except CMG4006 from La Asunta at 747 m) and has a wider geographical range from Coroico and La Asunta in the North to Cajuta and Quime in the South. Clade B is further resolved into a subclade (PP = 1.00) consisting of CMG4004 and CMG4215 and a main clade (PP = 0.88) with the remainder of the samples.

Alkaloid content in our tested samples was not significantly different among the three different trunk heights (Figure 2d) (Kruskal-Wallis test, quinine P = 0.40, quinidine P = 0.46, cinchonine P = 0.32, and cinchonidine P = 0.24) and consequently the average of the three values for each of the four major alkaloids was used in our subsequent analyses.

Total major alkaloid content varied significantly (Wilcoxon test = 99, P < 0.001) between the high alkaloid containing clade A (ranging from 13.6 to 41.3 mg g⁻¹, with a mean of 21.7 mg g⁻¹), and the lower alkaloid containing clade B (ranging from 0.3 to 12.5 mg g⁻¹, with a mean of 4.7 mg g⁻¹). Within clade B, one sample CMG4004 had a higher content that was in the range of samples from clade A, with a total major alkaloid content of 36.2 mg g⁻¹. Total major alkaloid variation across all samples varied from 0.3 mg g⁻¹ (CMG4218, <0.1% of bark) to 41.3 mg g⁻¹ (CMG4214; 4.1%) with a mean of 11.6 mg g⁻¹ (1.2%), which is in the lower range of what has previously been reported in materials of natural origin (Nair, 2010).

Quinine was the most abundant alkaloid with a mean of 6.6 mg g⁻¹, ranging from 0.2 to 25.8 mg g⁻¹ (<0.1–2.6%). While there were only two samples without any detectable levels of quinine, almost 50% of all samples tested had less than 5.0 mg g⁻¹ of quinine (Table 4; Figure 4a).

Quinidine was the second most abundant of the alkaloids, accounting for an average of 2.3 mg g⁻¹ (0.2%) across samples. Quinidine was present in all samples and highly variable, although less so than quinine, ranging between 0.1 and 10.0 mg g⁻¹ (Figure 4b). Cinchonine accounted for an average of 2.2 mg g⁻¹ (<0.2%). Like quinine, its content varied greatly, with nine trees having no content at all and half of all sampled trees having less than 2.0 mg g⁻¹ (Figure 4c), with one further sample (CMG4214) having exceptionally high content of cinchonine (22.9 mg g⁻¹). Cinchonidine was the least abundant of the major alkaloids, comprising an average of just 0.5 mg g⁻¹ (<<0.1% of bark). Only 7 of the 22 trees sampled contained any detectable level of cinchonidine, ranging between 0.6 and 3.1 mg g⁻¹ (Table 4; Figure 4d).

Spearman's rank-correlation coefficients were calculated to investigate the relationship between the content of the individual alkaloids. Cinchonidine and quinine were highly autocorrelated (R² = 0.69, P < 0.001), while cinchonine and cinchonidine also correlated (R² = 0.49, P = 0.020). There were however no other significant correlations among the alkaloids, with quinidine not correlating with any of the other alkaloids.

All four major alkaloids found within C. calisaya showed substantial variability among samples (Table 5). A significant phylogenetic signal was found for quinine (K = 1.299, P = 0.046), cinchonidine (K = 1.719, P = 0.002), and total major alkaloid content (K = 0.960, P = 0.020), whilst there was no correlation between phylogeny and quinidine (K = 0.316, P = 0.299) or cinchonine (K = 0.372, P = 0.325) (Table 3). Given the contrast in alkaloid content visually identified between the highly supported clade A compared to the rest of the tree (clade B), t-tests were performed to test for significance between the two clades. Again, quinine (t = 3.343, P = 0.004), cinchonidine (t = −3.277, P = 0.017), and total major alkaloid content (t = −2.739, P = 0.021) were significantly higher within clade A, while quinidine (t = 0.148, P = 0.885) and cinchonine (t = −0.992, P = 0.356) were not.

Climatic PC1 was shown to account for 20.2% of the variation within alkaloid production while geographical distance (as PCNM1) accounted for 9.7% of the variation, with a cumulative 29.9% of alkaloid variation accounted for within the analysis (Table 6). The second geographical component (PCNM2) and soil properties (soil properties PC1 and PC2) and year of sampling did not significantly influence alkaloid production.

The climatic PC1 significantly correlated with the alkaloid composition and closer inspection of its loadings revealed that altitude and precipitation predominantly contributed to its formation. Climatic effects on the alkaloids were investigated directly on individual alkaloid and total major alkaloid content using Spearman's rank correlation coefficient. Despite autocorrelation of altitude and precipitation, only altitude significantly correlated with any of the alkaloid concentrations, with altitude correlating with quinine (R_s = −0.548, P < 0.001; Figure 5A), cinchonine (R_s = −0.473, P = 0.026; Figure 5C), and cinchonidine (R_s = −0.507, P = 0.016; Figure 5D), while no correlation with quinidine was found (R_s = −0.005, P = 0.982; Figure 5B). Temperature did not autocorrelate with either altitude or precipitation, nor did it correlate with any of the major alkaloids.

Quinine, cinchonidine, and total major alkaloid content significantly correlated with phylogeny, whilst altitude also correlated with quinine content. Therefore, PGLS was performed to account for phylogenetic autocorrelation, testing for an independent altitudinal effect on quinine, cinchonidine, and total major alkaloid content. However, significance for altitude was lost on quinine (P = 0.236) and remained non-significant for cinchonidine (P = 0.851) and total major alkaloid content (P = 0.135).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer FG-F declared a shared affiliation, though no other collaboration, with several of the authors CM, CB, CC, EH, SH, NR to the handling Editor, who ensured that the process nevertheless met the standards of a fair and objective review.