Phylogeny Predicts the Quantity of Antimalarial Alkaloids within the Iconic Yellow Cinchona Bark (Rubiaceae: Cinchona calisaya)

Considerable inter- and intraspecific variation with respect to the quantity and composition of plant natural products exists. The processes that drive this variation remain largely unknown. Understanding which factors determine chemical diversity has the potential to shed light on plant defenses against herbivores and diseases and accelerate drug discovery. For centuries, Cinchona alkaloids were the primary treatment of malaria. Using Cinchona calisaya as a model, we generated genetic profiles of leaf samples from four plastid (trnL-F, matK, rps16, and ndhF) and one nuclear (ITS) DNA regions from twenty-two C. calisaya stands sampled in the Yungas region of Bolivia. Climatic and soil parameters were characterized and bark samples were analyzed for content of the four major alkaloids using HPLC-UV to explore the utility of evolutionary history (phylogeny) in determining variation within species of these compounds under natural conditions. A significant phylogenetic signal was found for the content of two out of four major Cinchona alkaloids (quinine and cinchonidine) and their total content. Climatic parameters, primarily driven by changing altitude, predicted 20.2% of the overall alkaloid variation, and geographical separation accounted for a further 9.7%. A clade of high alkaloid producing trees was identified that spanned a narrow range of altitudes, from 1,100 to 1,350 m. However, climate expressed by altitude was not a significant driver when accounting for phylogeny, suggesting that the chemical diversity is primarily driven by phylogeny. Comparisons of the relative effects of both environmental and genetic variability in determining plant chemical diversity have scarcely been performed at the genotypic level. In this study we demonstrate there is an essential need to do so if the extensive genotypic variation in plant biochemistry is to be fully understood.


INTRODUCTION
Bark from Cinchona trees (Cinchona L., Rubiaceae) of the Andean mountain forests produce quinine alkaloids, which were the only effective treatment of malaria for more than four centuries (Honigsbaum, 2001;Kaufman and Ruveda, 2005). The medicinal value of Cinchona bark was first discovered in Loxa (now Loja, Ecuador) in the seventeenth century by Jesuit monks, and soon exports of different varieties of Cinchona pubescens Vahl (red bark) from South America to Europe were reaching half a million kilograms bark per year (Roersch van der Hoogte and Pieters, 2015). Import could not meet demand, and a quest began for the most productive source of Cinchona trees to establish plantations by the British, Dutch, and French empires. The Bolivian Cinchona calisaya Wedd. (yellow bark, Figure 1) proved to be the most productive species known to date (Greenwood, 1992;Nair, 2010).
C. calisaya is one of 23 species of trees in the genus Cinchona described to date, which produce varying amounts of alkaloids. The four major Cinchona alkaloids (quinine, quinidine, cinchonine, and cinchonidine) (Figure 2a) all possess antimalarial activity but have different pharmacological profiles (Taggart et al., 1948;Hill, 1963;Bruce-Chwatt, 1990). Since the first isolation of quinine in 1820 over 30 minor and less studied Cinchona alkaloids have been described from the genus (Kacprzak, 2013). Bark and roots are the main source of Cinchona alkaloids, whereas cinchophyllines are reported from leaves (Aerts et al., 1991). The site of production of the alkaloids has not been established. In addition, Cinchona type alkaloids have also been found in the related genera Ladenbergia Klotzsch and Remijia DC (Okunade et al., 2001;Ruiz-Mesia et al., 2005;Cosenza et al., 2013).
Cinchona bark and its alkaloids remained the most efficient treatment of malaria until the 1940s when chloroquine and other synthetic antimalarial compounds were developed (Newman et al., 2000;Kaufman and Ruveda, 2005). With the development of resistant malaria strains (Bruce-Chwatt, 1990) the quest for new antimalarial compounds is continuing, and the discovery of artemisinin from a Chinese herbal medicine based on Artemisia annua L. (Tu, 2011), was rewarded with the Nobel prize in medicine in 2015.
Quinine content does not only vary among species (Nair, 2010), but also among populations from different sites, complicating the identification of the most productive Cinchona barks (Townley, 1922;Holland, 1932). Natural variation in quinine content remains unexplained, as few studies have been conducted in natural habitats until now (Rusby, 1931;Hodge, 1946;Bruce-Chwatt, 1990).
Alkaloids are a large and varied family of nitrogen-containing natural products occurring widespread across several lineages of vascular plant species with a high degree of specificity of subtypes to plant lineages (Zulak et al., 2006). They are not essential for primary metabolism, but play a number of specialized roles within the plants. Through selective up or down regulation, alkaloids can vary from complete absence to very high concentrations among individuals of the same species (Moore et al., 2014). Alkaloid production is one of the primary mechanisms of plants response to environmental changes (Theis and Lerdau, 2003;Ramakrishna and Ravishankar, 2011), and as with other plant chemical defenses, have likely developed over different evolutionary timescales (Becerra et al., 2009). Evolutionary approaches have been successfully implemented in predicting plant phytochemical composition, accelerating the discovery, and exploration of plant-based medicines (Bohlin et al., 2010;Zhu et al., 2011;Rønsted et al., 2012). However, other studies have found inconsistency of specialized metabolite profiles at various taxonomic levels (e.g., Wink, 2003;Wink and Mohamed, 2003). While it is established that different plants produce different specialized metabolites, the underlying causes determining these differences remain unknown. Two major hypotheses have been proposed, as outlined below.
The "escape-and-radiate" hypothesis (ERH) predicts sequential cycles between plant and herbivore/pathogens, with plants increasing in chemical complexity over evolutionary time, and the evolution of novel traits that promote speciation are incremental and directional throughout the diversification process (Ehrlich and Raven, 1964;Berenbaum, 1983;Agrawal, 2007).
The "resource availability hypothesis" (RAH) predicts that plants will invest more in defense when the cost of tissue replacement is high (Janzen, 1974;Coley et al., 1985;Fine et al., 2006). Therefore, in harsher, nutrient poor environments, there will be a greater investment in physical and chemical defense mechanisms over evolutionary time. For example, sparsely vegetated harsh soil is associated with increase in environmental stresses from greater exposure to droughts and herbivores, and is also expected to increase investment in chemical defenses or other morphological or physiological adaptations such as crypsis, escape, mimicry, etc. (Strauss and Cacho, 2013;Cacho et al., 2015).
Contrary to the ERH, the RAH was developed using different species rather than a single species. Furthermore, there is little evidence of this process within species under natural conditions. Early experiments within plantations of C. calisaya in Java, outside the natural habitat (Winters et al., 1947;Loustalot and Winters, 1948), found alkaloid content to increase with increasing soil nitrogen, while low levels of soil moisture caused a decrease in alkaloids, a clear contrast to the RAH (Arechavaleta et al., 1992;Malinowski et al., 1998).
Variation within species of plant chemical defenses has been the source of extensive investigations (Andrew et al., 2007;Barton and Koricheva, 2010;Kim et al., 2011;Holeski et al., 2012;Weldegergis et al., 2015). However, these studies have been mostly carried out in the laboratory under controlled conditions, reducing the complexity of interactions and variables that exist in native natural habitats (Moore et al., 2014). Additionally, the effect of genotypic variation within species on chemical profiles has not been studied extensively (Bidart-Bouzat and Kliebenstein, 2008). When these have been performed, genotypes have often been considered populations from separate geographical locations, leaving the genetic variation among individuals untested. Whether plant genotypes vary in their chemical defenses in Cinchona remains untested, although genotypic variation could potentially explain the substantial variation in alkaloids found within C. calisaya (Townley, 1922;Holland, 1932).
In addition to the development over evolutionary time frames, chemical defenses against pathogens and herbivores may be induced over a matter of hours or even minutes in response to external stimuli such as herbivory (Boege and Marquis, 2005;Barton and Koricheva, 2010) or over longer timespans in response to changing abiotic conditions (e.g., soil properties, climate, altitude) or biotic interactions (other plants, microorganisms) (Karban and Baldwin, 2007;Moore et al., 2014). The exact role of Cinchona alkaloid production in C. calisaya is unknown, but has been shown in laboratory experiments to inhibit larval growth (Green et al., 2002). It is thus hypothesized that quinine and the other Cinchona alkaloids are up-regulated in the presence of insects to inhibit their potential harmful activity on the trees.
The goal of this study is to explore the importance of genetic and environmental variation (such as climate, soil properties, and geographical distance) in determining variation of plant chemical defenses under natural conditions in C. calisaya. We also assess the applicability of the RAH in explaining resource allocation to chemical defense within a single species.

Collection of Cinchona calisaya Samples in Bolivia
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).

Detection of Alkaloids in the Bark
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 FIGURE 3 | Bayesian phylogram (50% majority-rule consensus) of the 22 specimens of Cinchona calisaya based on plastid (trnL-F, matK, rps16, and ndhF) and nuclear (ITS) sequences. The circle coloration represents the concentration of the total Cinchona alkaloids found in each sample on the tree and on the map ( Table 4). The reference map shows the distribution of the sampled species.
water (10:90 v/v) and mobile phase B: 100% methanol at a flow of 0.2 mL min −1 . 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).

DNA Sequencing and Phylogenetic Inference
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 2 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 partitionspecific 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 burnin phase (in this case the first 10,000 sampled trees).

Phylogenetic Signal in Alkaloid Production
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.

Environmental Variables
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. Alkaloid content are averages of samples at three heights per tree (0.7, 1.3, and 2.0 m). Alkaloid content in our tested samples was not significantly different among the three different trunk heights (Figures 2b,d).
Frontiers in Plant Science | www.frontiersin.org 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 2 formed by combustion in a muffle furnace (Nabertherm, Germany) in a CO 2 free environment (LECO Tru Mac N, Michigan). Total nitrogen was determined by LECO Tru Mac N as ammoniacal nitrogen NH 4 -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.

Multivariate Analyses of Explanatory Parameters
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.

Disentangling Phylogenetic and Environmental Variation on Alkaloids
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.

Variation in Climatic and Soil Parameters
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 −1 , 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 −1 ) were found for P. Both K and Mg varied considerably, K between 3.9 and 26.0 mg kg −1 , and Mg between 2.0 and 88.0 mg kg −1 among the sites ( Table 4).

Phylogenetic Structure of Cinchona calisaya Samples
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 −1 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.

Variation in Composition and Quantity of Major Antimalarial Alkaloids
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 −1 , with a mean of 21.7 mg g −1 ), and the lower alkaloid containing clade B (ranging from 0.3 to 12.5 mg g −1 , with a mean of 4.7 mg g −1 ). 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 −1 . Total major alkaloid variation across all samples varied from 0.3 mg g −1 (CMG4218, <0.1% of bark) to 41.3 mg g −1 (CMG4214; 4.1%) with a mean of 11.6 mg g −1 (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 −1 , ranging from 0.2 to 25.8 mg g −1 (<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 −1 of quinine (Table 4; Figure 4a).
Quinidine was the second most abundant of the alkaloids, accounting for an average of 2.3 mg g −1 (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 −1 (Figure 4b). Cinchonine accounted for an average of 2.2 mg g −1 (<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 −1 (Figure 4c), with one further sample (CMG4214) having exceptionally high content of cinchonine (22.9 mg g −1 ). Cinchonidine was the least abundant of the major alkaloids, comprising an average of just 0.5 mg g −1 (<<0.1% of  Table 4. bark). Only 7 of the 22 trees sampled contained any detectable level of cinchonidine, ranging between 0.6 and 3.1 mg g −1 (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 2 = 0.69, P < 0.001), while cinchonine and cinchonidine also correlated (R 2 = 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.

Multivariate Analyses of Explanatory Parameters
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.

Disentangling Phylogenetic and Environmental Effects on Alkaloid Production
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 Bold indicates significance using Spearman's rank-order correlation (P < 0.05). PC1 is the principal component for the climatic parameters, PCNM1 and PCNM2 are principal coordinates of neighbors for geographical distances, and Soil properties PC1 and 2 were the first 2 principal components for the soil properties. for cinchonidine (P = 0.851) and total major alkaloid content (P = 0.135).

DISCUSSION
Historically, quinine and quinine-like alkaloid concentrations are known to vary considerably within C. calisaya (Delondre and Bouchardat, 1854), however what is driving this variation is not known. This work represents one of the first attempts to simultaneously explore the relative effects of genotype and environmental variation in determining plant secondary metabolite production under natural conditions. Our results show that genotypic variation within C. calisaya is the primary driver of differences in alkaloid content, rather than environmental variation. A significant correlation was found between phylogeny and quinine, cinchonidine, and total major alkaloid production, whilst there was no significance for quinidine and cinchonine. Within the phylogeny, two highly supported clades were identified (clade A and clade B). Clade A was a small mostly geographically clustered group that contained all but one of the higher alkaloid producing samples, whilst clade B was comprised of a number of further subclades. Clades A and B differed significantly in alkaloid content and it is likely differences between these two drove the correlation between phylogeny and alkaloids. Additionally, the link between phylogeny and alkaloid content may have been confounded within this work. Both the high producing sample in clade B found in Chulumani (CMG4004) and the low producing sample in clade B found in Coroico (CMG4215) were grown around settlements and were possibly cultivated plants that could have been transplanted recently from other areas, which would explain the inconsistency between geography and phylogeny.
Our results also suggested that clade A was geographically clustered, centered between Chulumani and Coripata in Bolivia. Plants are more likely to be related the closer they are geographically, therefore it is unsurprising that geographical separation correlated with alkaloid composition given the significant phylogenetic effect (Castells et al., 2005;Egydio et al., 2013). They are also more likely to share similar environmental conditions and samples in clade A were also found at lower altitudes than the other samples, ranging between 1,100 and 1,350 m. With the exception of a single sample (CMG4006), which was found at 747 m, all other samples from clade B were found above 1,700 m and produced less of the major alkaloids than the other samples in clade B. Climatic factors are thought to be amongst the most important factors determining alkaloid composition in plants (Lobo et al., 2010;Moore et al., 2014), and the negative correlation between altitude and chemical defense compounds which we found for Cinchona spp. has been observed before in for example the case of Artemisia (Fluck, 1963). However, to our knowledge a comparison between genotypic variation and environmental effects on alkaloid chemistry in the field has not been previously investigated, and in this study environmental effects were no longer significant once phylogenetic effects were taken into account. Results here therefore suggest the development of a locally adapted, highly productive genotype of C. calisaya in the mountain forests of western Bolivia over evolutionary rather than ecological time scales. This could therefore be compatible with the ERH (Ehrlich and Raven, 1964;Berenbaum, 1983;Agrawal, 2007), in which genotypes eventually become separate species with differing chemotypes.
Changing altitude will also affect a number of environmental properties (Del Grosso et al., 2008;Sundqvist et al., 2013), with higher rates of photosynthesis associated with lower altitudes that are generally warmer (Berry and Bjorkman, 1980;Korner and Diemer, 1987). Meanwhile, no significant effects of soil properties were observed on the four major alkaloids under natural conditions. This is in contrast to experimental manipulations, which suggest a link between alkaloid production and soil phosphorus, nitrogen and pH (Koblitz et al., 1983;Harkes et al., 1985;Arechavaleta et al., 1992;Malinowski et al., 1998), although soil properties varied under natural conditions by orders of magnitude less than the conditions established experimentally. With no significant effects of nutrient availability and potentially enhanced rates of photosynthesis of samples of clade A found in lower altitudes, our results suggest that the RAH may not apply within species despite evidence of a phylogenetic effect.
In contrast to previous studies (Nair, 2010), no significant variation was found between barks sampled at heights of 0.7, 1.3, and 2.0 m. Generally, the alkaloids are thought to be in lower content in the twigs, higher in trunk bark, and maximum in the root bark (Nair, 2010) with peak quinine content found in the collar portion (30-45 cm in length, near the base) (Hodge, 1946;Nair, 2010). However, we found limited evidence for variation in any of the four major alkaloids across the bark of the trunk. A significant proportion of alkaloid variation remained unexplained, particularly within quinidine and cinchonine.
Although temporal and biotic conditions were not studied in this work, future studies may greatly benefit from this, with sampling under natural conditions assisting in discovering novel mechanisms regulating plant secondary metabolites. For example, microbial communities in the bark have shown ability to influence production of Cinchona alkaloids (Koblitz et al., 1983;Maehara et al., 2010Maehara et al., , 2011 but only in vitro. Furthermore, quinine and quinine like alkaloids have demonstrated activity against larvae (Green et al., 2002) in vitro, but no study to date has investigated the role of insect herbivory on C. calisaya alkaloid production using whole plants. Interestingly, insect activity and composition are known to be variable across altitudinal gradients, which peaks over the intermediate altitudes that clade A occurred, decreasing at the higher altitudes clade B was found (Rahbek, 1995(Rahbek, , 2005Kessler et al., 2011;Larsen et al., 2011;Pellissier et al., 2012). Further studies understanding insect biodiversity and composition could therefore be a fruitful area of research in understanding alkaloid regulation (Beccera, 2015). It should however be noted that studies under natural conditions can only imply causation and generate hypotheses, and controlled environment studies should be employed in combination with environmental sampling to maximize our understanding of plant secondary metabolite production. In conclusion, we hypothesize that phylogenetic variation is the main driver of changes in alkaloid composition. Growing plants from both clade A and B in controlled environments may help confirm this hypothesis.

AUTHOR CONTRIBUTIONS
NR and CM conceived the ideas and supervised the research. SH, CC, and EH conducted the chemical analysis. CM conducted the fieldwork. CM assembled and analyzed the data with contributions from CB. CB designed and conducted the statistical analyses. CM wrote the manuscript with contributions from CB, NR and AA. All authors read and improved the manuscript.