Root Endophytes of Coffee (Coffea arabica): Variation Across Climatic Gradients and Relationships with Functional Traits

It has been suggested that both plants and animals should be seen as holobionts, composite meta-organisms comprised of hosts and their microbiomes that experience environmental challenges as a unit (Agler et al. 2016; Vandenkoornhuyse et al. 2015; Zilber-Rosenberg and Rosenberg 2008). The diversity of both fungi and bacteria that reside within plant tissues (endophytes) or in close proximity to them in the soil (rhizosphere communities) is proving to be much greater than previously appreciated (Bulgarelli et al. 2013; Vandenkoornhuyse et al. 2002). The importance of mycorrhizal fungi for enhancing water and nutrient acquisition by host plants is well known (Bonfante and Anca 2009; Koltai and Kapulnik 2010) and they are recognized as playing a key role in plant and ecosystem resistance and resilience to environmental change (Furze et al. 2017). Bacterial endophytes are equally influential, playing key roles in plants including the production and/or modulation of plant growth regulators (Santoyo et al. 2016), the provision of nutrients via non-rhizobial N₂ fixation (Asis and Adachi 2004; Pham et al. 2017; Ruppel et al. 1992) and phosphorus solubilization (Padder et al. 2016; Walia et al. 2017), the provision of plant protection from pathogens via production of antibiotic compounds (Compant et al. 2005), and enhanced plant salt tolerance (Ali et al. 2014).

A number of studies have contributed insights into the determinants of plant microbiome community structures. For example, host tree taxonomic identity explained over 50% of microbial community composition in the phyllosphere microbial communities of coexisting tropical tree species (Kembel and Mueller 2014; Kembel et al. 2014). Similarly, Fitzpatrick et al. (2018) found root endosphere community similarities were positively correlated to phylogenetic relatedness between plant species in angiosperm for root bacterial endophyte communities. Other large amounts of variance are determined by plant tissue type or compartment (i.e., root versus stem or leaf, Coleman-Derr et al. 2016), growing season (e.g., Shen and Fulthorpe 2015), host soil conditions or growth media (Bulgarelli et al. 2012, 2013; Gottel et al. 2011; Lundberg et al. 2012; Ringelberg et al. 2012), or in response to temperature and precipitation change (Castro et al. 2010; Coleman-Derr et al. 2016). What remains not fully understood is the extent to which plants have a core microbiome of ecophysiological importance or just microbial communities that are stochastic assemblages drawn from nearby soils and the atmosphere. The presence of a core microbiome would imply that these communities are structured by plant-specific selective factors, and would not only drive future studies on their main interactions, but potentially be optimized for greater ecosystem function (Toju et al. 2018).

The definition of core microbial communities in both crop and noncrop plants has been identified as a key research priority for our increasingly crowded planet (Busby et al. 2017). Coffee (Coffea spp.)—the focus of our study here—is among the world’s most widespread crops, covering nearly 11 million ha of agricultural land (FAO 2018). Coffee microbial communities have been examined for their ability to protect plants from pathogens (e.g., the devastating coffee leaf rust, Hemileia vastatrix Berk. & Br. [Shiomi et al. 2006; Silva et al. 2012]) and for their influence on coffee bean production and quality (Oliveira et al. 2013; Vaughan et al. 2015). Studies focused on the relationship between microbial associates and coffee management systems are more recent. For instance, variation in arbuscular mycorrhizal fungi and bacterial endophyte communities across gradients of management intensification has been evaluated, with indications that root arbuscular mycorrhizal fungi decline dramatically in monocultures relative to more diversified production systems (De Beenhouwer et al. 2015b; Munroe et al. 2015). Caldwell et al. (2015) used pyrosequencing of bacterial sequences to compare coffee rhizosphere to control soils at organic, intensive, and transitional coffee plantations. Each site proved to have a unique community structure, particularly with respect to nitrogen fixing genera.

Also remaining near-completely unexplored is how changes in endosphere microbial diversity within a host species across environmental gradients correlates with plant functioning. Specifically, the plant functional ecology literature suggests that plant response to environmental change, and impacts on rates of ecosystem function, are mechanistically governed by a relatively small number of plant functional traits (Diaz et al. 2016; Reich 2014; Reich et al. 1997, 1999; Violle et al. 2007; Westoby 1998, 2002). Relationships between functional traits and belowground microbial diversity remain limited; indeed, elucidating such relationships represent one of the remaining frontiers in trait-based ecology (Laliberte 2017). Arguably the most widely studied plant traits are those comprising the leaf economics spectrum: six covarying leaf traits including leaf mass per area, maximum photosynthetic (A) and dark respiration rates (R_d), leaf N and P concentrations, and leaf lifespan (Reich et al. 1997, 1999; Wright et al. 2004, 2005). Multiple studies have shown that inter- and intraspecific variation in leaf economics spectrum traits are key determinants of ecosystem functions including phyllosphere microbial community composition (Kembel and Mueller 2014; Kembel et al. 2014) and plant susceptibility to pests or pathogens (Poorter et al. 2004). Authors have also demonstrated that variation in these same leaf economics spectrum traits within crop genotypes is critically important in governing rates of agroecosystem functions such as yield and nutrient cycling (reviewed by Martin and Isaac 2015, 2018).

Analogous to the leaf economics spectrum, researchers have demonstrated that plant roots differentiate from one another along a root economics spectrum that is linked with belowground resource acquisition strategies among and within plants (Hajek et al. 2013; Isaac et al. 2017; Ostonen et al. 2007). While a definitive account of the traits forming the root economics spectrum continues to evolve, the root economics spectrum is generally hypothesized to consist of specific root length, specific root area, specific root tip density, and root diameter (Laliberte 2017; Ma et al. 2018; Prieto et al. 2015; Roumet et al. 2016; Weemstra et al. 2016). In turn, root economics spectrum traits are key predictors of water and nutrient acquisition (Prieto et al. 2015), rates of opportunistic root proliferation (Eissenstat 1991), the C economy of roots vis-à-vis root respiration and decomposability (Roumet et al. 2016), and plant responses to soil fertility and moisture (Isaac et al. 2017). Friesen et al. (2011) reviewed hypotheses on how certain root economics spectrum traits may be mediated by microbial diversity. While their analysis presents expectations of microbial mediation of some root economics spectrum traits, they note that evidence supporting these relationships is generally weak, limited to a small number of host herbaceous species, and pertains only to a very small number of specific host−microbe relationships (Table 1 in Friesen et al. 2011).

Our study was designed to explore intraspecific variation in the root endosphere microbiome, and better understand what explains this variation. Using the economically important coffee tree crop along an environmental gradient and under contrasting management conditions, our study sought to address the following questions. (i) Is there a core coffee root microbiome? If so, (ii) do environmental (edaphic and climatic) and management (shade grown and monoculture) conditions explain variation in this microbiome? And (iii) do plant functional traits and suites of traits covary with the root microbiome?

MATERIALS AND METHODS

Study site and design.

Our study assessed bacterial and fungal root endophyte communities in coffee across a wide climatic gradient within the coffee growing regions of Costa Rica and Nicaragua (Fig. 1A). Within these countries, four different study sites were selected in order to represent the broad range of coffee growing conditions in terms of temperature and precipitation (Fig. 1B and C). These sites included the following: (i) a hot and wet site (The Centro Agronomico Tropical de Investigacion y Ensenanza [CATIE]; mean annual temperatures [MAT] 23.4°C, and mean annual rainfall [MAR] ∼3,200 mm); (ii) a cool and wet site (Aquiares; MAT 19.5°C, MAR ∼3,000 mm); (iii) a cool and dry site (Llano Bonito; MAT 18.7°C, MAR ∼1,500 mm); and (iv) a hot and dry site (Masatepe; MAT 24°C, MAR ∼1,400 to 1,470 mm) (Fig. 1A to C). At these sites two coffee cultivars are grown, namely C. arabica ‘Caturra’ at CATIE, Aquiares, and Llano Bonito, and C. arabica ‘Pacas’ at Masatepe, which is a genotype long recognized as a mutation of the same genetic strain as C. arabica ‘Caturra’ (Bertrand et al. 1999). The optimal range for coffee performance is MAT of 18 to 21°C (Lin 2007) and MAR of 2,200 mm (Vaast et al. 2006). Therefore, CATIE and Masatepe represent higher than optimal temperatures, while Llano Bonito and Masatepe sites receive suboptimal rainfall.

Download as PowerPoint

We made use of the two distinct management arrangements implemented at all sites including (i) full sun monoculture, and (ii) shade coffee agroforestry where coffee is intercropped with N₂-fixing leguminous tree species, namely, Erythrina poeppigiana (Walp.) O.F. Cook (Fabaceae) (at CATIE, Aquiares, and Llano Bonito) and Inga laurina (Sw.) Willd. (Fabaceae) (at Masatepe). Within each site, two separate 25-m² research blocks were identified within both management conditions, situated a minimum of 20 m apart. Within each block, we selected four coffee plants for evaluation of belowground microbial community composition under the following design: n = 64 coffee plants situated evenly across n = 16 sample blocks, which were located within n = 2 different management systems across n = 4 sites.

Coffee plant selection and environmental parameters.

In order to account for potentially confounding effects of plant ontogeny and seasonality on microbial communities, sampling took place during a 4-week period in April 2014, and all 64 individual coffee plants selected for analyses were reproductively mature stump-pruned stems, measuring 140 to 235 cm in height with resprout basal diameters between 14.4 to 34.6 cm. Plants in the agroforestry systems were situated within 0.5 to 9.6 m from the nearest shade tree, while those in the monoculture systems were a minimum of 20 m from the nearest shade tree to avoid interacting effects such as N transfer (Martin et al. 2017). For each individual plant, we also measured a suite of environmental variables. Field and lab methods associated with collection of environmental data are described in detail in our previous analyses of coffee functional traits (Isaac et al. 2017; Martin et al. 2017), but in short we assessed canopy openness (%), soil moisture content (%), soil pH, soil C and N (%), and plant available P (mg kg⁻¹) for each individual plant.

Collection of leaf and root functional traits.

Coffee leaf samples were collected following protocols detailed by Martin et al. (2017) that were based on well-established trait collection protocols (Perez-Harguindeguy et al. 2013). The following leaf traits were measured on six leaves per the plant and plant-level means used for our analysis here: mass-based photosynthesis (A_mass, μmol CO₂ m⁻² s⁻¹), water use efficiency (mmol CO₂ mol H₂O⁻¹), leaf density (g cm⁻³), leaf mass per area (g m⁻¹), leaf N (% mass), and leaf area (cm²). Root functional traits were derived from root samples collected following protocols detailed by Isaac et al. (2017), which were informed by established root trait collection protocols (McCormack et al. 2015; Perez-Harguindeguy et al. 2013). The following root traits were calculated for each plant: root diameter (mm), specific root length (m g⁻¹), specific root area (m² kg⁻¹), specific root tip density (number of tips g⁻¹), root length density (cm cm⁻³ of soil), root N (% mass), and root C:N. Root traits were measured based on one complete intact lateral root that was excised using a sterilized scalpel, immediately placed in sealed polyethylene bags, and stored at −20°C within 2 h of collection until processing at the University of Toronto Scarborough. The same roots used for root trait determinations were subsequently used for assessments of microbial communities. Root tips (2 cm section from root ends) were collected, stored, and analyzed separately from proximal roots (2 cm sections collected from lateral roots in the top 20 cm).

Processing of roots.

In the lab, roots were first superficially cleaned of soil by rinsing in tap water, and subsequently carefully washed in 0.1% Triton X-100 solution, rinsed in distilled water, and surface sterilized through successive washes in 70% ethanol and 1.5% bleach. To assess the effectiveness of the surface sterilization, the final sterile water wash was spread plated, and the sterilized surfaces were imprinted onto agar plates.

From each sample, between 175 and 200 mg of wet roots were used for DNA extraction using FastDNA Spin Kits (MP Biomedicals, Lachine, QC, Canada). The extracted DNA was quantified using fluorescence estimates from gel electrophoresis, and from UV absorbance as determined on a NanoDrop Spectrophotometer. Fast DNA Spin Kits were extracted without samples in order to detect contaminants. The resulting aliquots were subjected to amplification and terminal restriction fragment length polymorphism (TRFLP) analyses. These blank samples did not produce any detectable amplicon fragments.

TRFLP on all samples.

In order to assess the variability of root microbial communities between trees in a cost-effective manner, we analyzed all of the root samples using TRFLP. For bacterial community analysis via TRFLP, rDNA fragments were amplified from samples using primers 27f-FAM and 1492r-HEX in conjunction with blocker primers used to minimize amplification of plant cell organelle rDNA fragments (Shen and Fulthorpe 2015). Amplicons were digested using MspI and sent to the Genomics Facility of the Advanced Analysis Centre, University of Guelph, Canada, for fragment analysis on an Applied BioSystems 3730 DNA Analyzer. For fungal community analyses, ribosomal interspacer fragments were amplified using internal transcribed spacer (ITS)1F 5′-CTT GGT CAT TTA GAG GAA GTA A-3′ forward labeled with phosphoramidite 6-FAM and ITS4 5′-TCC TCC GCT TAT TGA TAT GC-3′ reverse labeled with Cal Red. The amplicons were digested in BfaI prior to fragment analysis. For both bacterial and fungal datasets, terminal fragments less than 50 bp were deleted and any terminal fragments (phylotypes) totaling less than 5% across all samples were not included in subsequent analyses.

Sample pooling, predigestion, and sequencing.

In order to obtain identification of the taxa present in our samples in an affordable way, we combined the tree root tip DNA samples derived from individual trees into eight pools, one for each site/management group. As we also had proximal root DNA from the Aquiares and CATIE sites, we made an addition four pools from these (a sun and shade pool each at each site). To generate the pools, the same mass of DNA was taken from each subsample and mixed to ensure no one coffee tree was replicated or overrepresented. For bacterial analyses, a subsample of each pool was digested with PvuII and MscI prior to shipment to MR DNA labs (Shallowater, TX, U.S.A.) for amplification and Illumina sequencing. Digestion with these enzymes limits the amount of amplification of ribosomal genes from the plant organelles (Shen and Fulthorpe 2015). This step was not required for the analysis of fungal ITS sequences.

At the MR DNA lab, samples were amplified using ITS1F and ITS4 for fungal analysis and 27f and 806r for bacterial analyses. Sequencing was performed on a MiSeq following manufacturer’s guidelines. Sequence data were then joined, depleted of barcodes, and removed from database if they were <150 bp in size or had ambiguous base calls; the remaining sequences were denoised, and chimeras were removed. Operational taxonomic units (OTUs) were defined by clustering at 3% divergence (97% similarity). Final OTUs were taxonomically classified using BLASTn against a curated database derived from RDPII (rdp.cme.msu.edu) and NCBI (https://www.ncbi.nlm.nih.gov/).

Statistical analyses of sequence data–site level community compositions.

Sequence data on all bacterial and fungal OTU abundances were reduced as follows. Data were first converted to percentage of total reads for each pool. Where we had both root tip and proximal root data, they were highly correlated to each other (r > 0.95) so the data were averaged for these two compartments. The most numerically abundant strains were ordered for each site and management treatment and compared. Bacterial OTUs that were present at all four sites and averaged 0.5% or more in abundance were retained as core. Fungal OTUs varied more between sites, so those present at three sites or more and found within the top 30 ranked abundances were retained as core, as were three of the most abundant arbuscular mycorrhizal fungi (AMF) species. This data reduction procedure resulted in 26 bacterial genera and 31 fungal genera being retained for a network analysis.

Specifically, we assessed relationships among these bacterial and fungal genera abundances in coffee across all eight site-by-management conditions using Pearson correlations performed on arcsine-transformed relative abundances of each genera in our core dataset (represented by a given OTU classified to genus, and where n = 8 for all pairwise correlations). To account for possible spurious correlations occurring as a result of the large number of comparisons (i.e., 3,192 tests in total), the significance of all correlation coefficients (Pearson’s) was evaluated against a type 1 error rate of 0.01. The correlation network of the core microbiome was then determined using network analysis visualization techniques implemented with the ‘qgraph’ function in the ‘qgraph’ R package (Epskamp et al. 2012).

Statistical analyses–Intrasite variation and differences between sites and management treatments.

All data on phylotypes (terminal fragments) obtained from the TRFLP analyses, i.e., our tree level data, were used to evaluate if bacterial and fungal communities (assessed separately) differed significantly as a function of site, shade management conditions, a site-by-shade interaction term. This was done using nonmetric multidimensional scaling (NMDS) based on Bray-Curtis dissimilarity values among all samples, which were subject to a Wisconsin double standardization transformation prior to analysis. All NMDS analyses were implemented using the ‘metaMDS’ function in the ‘vegan’ R package (Oksanen et al. 2016). We then used a permutational analysis of variance (PERMANOVA) implemented using the ‘adonis’ function in the vegan R package, to evaluate if microbial community dissimilarities were significantly related to site, shade treatment, a site-by-shade treatment interaction term. In short, these analyses indicated that bacterial and fungal communities differed significantly as a function of site and site-by-shade interaction, which then informed our next analysis.

We used redundancy analyses (RDA) to evaluate if differences in bacterial and fungal communities among samples were significantly predicted by continuous environmental variables and/or leaf and root functional traits. For both bacterial and fungal communities, we performed separate RDAs, which included the following predictors: (i) seven environmental variables (canopy openness, soil P, soil pH, soil N, soil C, and soil moisture); (ii) seven root functional traits (root diameter, specific root length, specific root area, specific root tip density, root length density, root N, and root C:N); and (iii) six leaf functional traits (A_mass, water use efficiency, leaf density, leaf mass per area, leaf N, and leaf area). These three sets of environmental variables and traits were included as variables in the constrained RDA, whereas site and shade treatment were included as actors (as per the results of our NMDS analyses). All RDAs were performed based on Hellinger-transformed phylotype abundance data (as suggested by Legendre and Gallagher [2001]), and significance of each variable was assessed using a PERMANOVA implemented using the ‘anova.cca’ function in the ‘vegan’ R package (with 9,999 permutations used).

RESULTS

The core microbiome–Microbial genera present in roots and their correlations.

Illumina sequencing of ribosomal amplicons provided 17,437 to 46,321 ITS (fungal) sequence reads representing 54 to 112 fungal genera, and 15,110 to 74,088 16S (bacterial) reads representing 218 to 319 bacterial genera from each of the DNA pools (Supplementary Table S1). As the root tip and proximal root data were highly correlated for each site/management pair, we averaged the data from those root sources for the CATIE and Aquiares sites. Sequence data are archived in the sequence read archive (SRA) under BioProject accession number PRJNA556713.

Across all eight site-by-management system study locations, the highest mean relative abundances in the root bacterial root communities were seen for Pantoea (17.2%), Enterobacter (12.2%), and Burkholderia (7.6%). On average, the most abundant fungal genera were Cladosporium (40%), Penicillium (12%), Exidiopsis (6.6%), and Trechispora (5.4%) and Mycena (3.9%) (Table 1, Fig. 2).

TABLE 1 Relative abundances of core endophytic bacterial and fungal genera in coffee (Coffea arabica), across four growing conditions and two management treatments^a

		Cool, wet (Aquiares)		Hot, wet (CATIE)		Cool, dry (Llano Bonito)		Hot, dry (Masatepe)
Endophyte group	Site management	FS	SH	FS	SH	FS	SH	FS	SH	All sites
Bacteria	Acidipila	0.37	0.85	0.63	0.18	0.44	0.04	2.30	1.87	0.84
	Acidobacterium	1.33	2.33	1.45	0.08	1.75	0.27	2.14	3.97	1.67
	Actinoallomurus	0.27	0.56	0.31	0.07	0.55	0.09	0.72	2.73	0.66
	Amycolatopsis	1.46	2.97	2.03	0.62	3.99	0.66	6.17	5.33	2.90
	Blastopirellula	1.40	3.20	0.50	0.06	1.03	0.26	1.13	1.39	1.12
	Bradyrhizobium	1.35	1.86	1.82	0.38	7.95	1.62	1.54	2.91	2.43
	Burkholderia*	3.38	3.85	2.79	0.52	18.59	6.33	13.30	12.38	7.64
	Edaphobacter	3.70	5.79	2.10	0.20	0.75	0.19	0.99	4.29	2.25
	Enterobacter	21.83	11.18	4.66	10.70	0.37	48.39	0.42	0.11	12.21
	Erwinia	4.26	2.90	0.58	0.21	0.27	14.99	0.37	0.15	2.97
	Frankia	1.49	2.35	2.41	0.70	3.16	0.40	3.87	2.67	2.13
	Gemmata	1.28	1.93	0.43	0.09	3.25	0.53	0.64	4.45	1.57
	Granulicella	0.39	0.53	0.15	0.02	0.07	0.02	4.22	0.64	0.75
	Humibacter	1.94	5.41	2.40	1.36	0.09	0.03	0.56	1.64	1.68
	Ktedonobacter*	0.46	0.85	1.15	0.29	9.79	1.80	13.42	6.27	4.25
	Kutzneria	0.04	0.06	0.07	0.03	0.63	0.24	3.89	0.11	0.63
	Mycobacterium	1.20	2.62	1.17	0.49	0.94	0.26	1.72	4.27	1.58
	Nitrosovibrio	0.58	1.50	0.53	0.13	0.77	0.25	0.68	0.81	0.65
	Nocardia	0.39	0.77	0.72	0.13	1.55	0.38	8.62	3.38	1.99
	Pantoea*	19.85	5.33	43.54	53.91	1.21	12.05	1.62	0.39	17.24
	Pseudomonas*	1.17	0.85	5.57	17.80	0.23	2.61	0.29	0.10	3.58
	Ralstonia	1.09	2.26	0.23	0.03	0.55	0.21	0.87	2.26	0.94
	Saccharibacter	0.90	1.28	1.92	0.70	0.66	0.08	0.84	0.34	0.84
	Solirubrobacter	1.82	1.72	1.79	0.64	0.71	0.28	0.36	0.28	0.95
	Steroidobacter	1.38	3.07	1.99	0.22	2.87	0.49	4.05	7.40	2.68
	Thermosporothrix	0.05	0.06	0.13	0.04	0.15	0.08	4.10	0.21	0.60
Fungi	Acaulosporab	0.10	0.01	0.26	0.01	0.00	0.00	0.00	0.00	0.05
	Acremonium	1.28	2.30	0.14	0.37	0.04	0.12	0.15	0.01	0.55
	Archaeosporab	0.00	0.00	0.03	0.04	0.00	0.01	0.00	0.00	0.01
	Aspergillus	0.04	0.40	0.51	0.11	0.47	0.89	0.05	0.13	0.32
	Camarosporium	0.00	0.00	0.06	0.00	0.01	0.00	4.19	13.85	2.26
	Cladosporium**	49.45	50.67	47.66	42.20	81.29	47.44	0.17	0.17	39.88
	Clitocybe	0.00	0.04	0.00	5.00	0.00	0.00	0.00	0.00	0.63
	Codinaeopsis**	0.02	0.00	0.05	0.00	0.25	0.01	8.25	11.83	2.55
	Colletotrichum	1.63	0.84	3.56	0.26	0.01	0.00	0.01	0.01	0.79
	Conoplea	0.00	0.00	0.00	0.00	0.00	9.21	0.00	0.01	1.15
	Cryptococcus	1.91	2.10	1.90	0.71	0.60	3.08	1.72	7.37	2.42
	Dokmaia***	1.48	2.90	2.39	1.52	0.06	0.30	0.00	0.01	1.08
	Exidiopsis**	0.01	0.02	0.02	0.01	0.02	0.01	40.00	13.04	6.64
	Fusarium	15.64	4.14	2.32	3.63	0.43	0.96	0.37	2.29	3.72
	Fuscoporia	0.17	0.02	2.99	0.00	0.51	0.00	0.00	0.01	0.46
	Glomusb	0.01	0.01	0.03	0.05	0.01	0.05	0.20	0.00	0.05
	Isaria	0.16	0.91	0.03	0.02	0.00	0.00	0.00	0.00	0.14
	Knufia	0.01	0.00	0.00	0.00	0.00	0.00	4.90	8.49	1.68
	Marasmius	0.00	1.27	0.00	0.22	1.22	0.00	0.00	0.00	0.34
	Mycena	1.22	1.23	1.46	1.52	1.66	8.00	0.00	16.05	3.89
	Paecilomyces	0.63	0.63	0.43	1.01	0.06	0.12	0.00	0.09	0.37
	Paramicrosporidium	0.09	0.21	0.16	0.08	0.04	11.58	0.40	1.77	1.79
	Penicillium	14.72	6.44	30.48	29.34	7.08	5.98	1.02	0.73	11.97
	Phoma	4.38	15.14	1.06	2.93	1.45	1.24	0.01	0.03	3.28
	Plectosphaerella	0.16	0.58	0.02	0.38	0.01	0.05	0.00	0.00	0.15
	Stagonosporopsis	0.10	0.78	0.01	0.33	1.12	0.22	0.01	0.00	0.32
	Trechispora	0.00	0.06	0.01	1.52	0.01	0.01	30.89	10.91	5.43
	Trichoderma	2.48	0.10	0.07	0.93	0.16	0.41	0.06	0.25	0.56
	Trichosporon	0.00	0.09	0.00	0.01	0.10	0.00	0.00	0.01	0.03
	Verticillium	0.52	0.97	0.07	2.50	0.01	0.31	0.00	0.02	0.55

^a Climate data and locations for each site are presented in Figure 1A, and management treatments correspond to coffee growing in full sun (FS) monoculture and shaded (SH) agroforestry. Also shown are average relative abundances for each genera across all sites and management conditions. *, **, and *** indicate site differences are significant at 0.05 < P < 0.10, P < 0.05, and P < 0.001, respectively, according to two-way analysis of variances without replication. No differences between FS and SH at a given site were significant.

^b Genus is arbuscular-mycorrhizal.

View as image HTML

Download as PowerPoint

Only 172 of a possible 3,192 pairwise correlations performed among the relative abundances of the core 26 bacterial and 31 fungal species were significant (P ≤ 0.01; Supplementary Table S2). The network analysis (Fig. 2) shows three genera not correlated to any other: Glomus, Humibacter, and Aspergillus. There were two pairs uncorrelated to other taxa: Trichoderma and Fusarium, and Mycena and Cryptococcus. Three genera formed in a small correlated trio: Stagonosporopsis, Marasmius, and Trichosporon. The rest formed two loose networks connected through correlations to Ktedonobacter. Ktedonobacter was positively correlated to Amycolatopsis, Kutzneria, and Nocardia in a large, mostly bacterial, network on one side, and positively to Burkholderia, negatively to Paecilomyces and Dokmaia on the other. Abundant microbiome members Pantoea and Pseudomonas are positively correlated to each and to Penicillium and Archaespora. Enterobacter, Erwinia, Conoplea, and Paramicrosporidium form a tail connected to the majority via negative correlation to Frankia.

Overall, Nocardia, Acidipila, Codineopsis, and Exidiopsis exhibited the most positive correlations to other genera (9, 8, 7, and 7, respectively, versus an average of 2.6, Supplementary Table S3). Cladosporium was negatively correlated to five other fungal genera (Camarosporium, Codinaeopsis, Exidiopsis, Knufia, and Trechispora)—much higher than the average number of negative correlations (0.4).

The relative abundances of the 11 most abundant bacterial genera varied with climatic zone regardless of management practices (Fig. 2B). The wet sites had slightly higher abundances of Pantoea and Pseudomonas than the dry sites, which had higher levels of Burkholderia and Ktedonobacter (two-way ANOVA without replication, P < 10). Fungal relative abundance differences were greater with coffee roots at the hot and dry site, being dominated by Exidiopsis and Codinaeopsis at the dry sites, while the wetter sites were associated with much higher levels of Cladosporium and Dokmaia (Table 1).

Variation in coffee microbial communities across climate and management.

Tree level TRFLP analyses revealed a total of 211 fungal and 91 bacterial phylotypes in our data. NMDS analyses of these data (Supplementary Fig. S1) indicated that microbial community compositions differed significantly between sites in spite of high tree level intrasite variation (Adonis r² = 13.8 and 11.5 for bacterial and fungal communities, P < 0.001 in both cases, Supplementary Table S4). While there was also a significant site-by-management interaction effect for both bacteria and fungi (Adonis r² = 7.9 and 2.2, respectively, P ≤ 0.022), microbial communities did not differ systematically between sun and shade management systems (Adonis test r² = 0.2 and 1.3 for fungi and bacteria, respectively, P ≥ 0.151).

Variation in bacterial communities across environmental gradients and with functional traits.

RDA analyses indicated that bacterial community composition varied significantly as a function of canopy openness, soil pH, and soil C concentrations (P ≤ 0.02 in all three cases) (Table 2). Specifically, variation in bacterial community composition along the first RDA axis was related to differences in soil pH among sites, with the two dry sites expressing high soil pH values compared with the two wet sites (Fig. 3A). Bacterial community variation along the first RDA axis was also linked with site differences in soil carbon, with these values being highest in the cool and wet site (Fig. 3A). Bacterial composition along the second RDA axis was related to site differences in canopy openness, primarily driven by high openness values in the hot and dry site (Fig. 3A). Sites also differed in terms of the coffee root traits values being expressed in patterns generally consistent with hypotheses of a root economics spectrum: plants in the dry sites expressing roots with more resource conserving traits versus more resource acquisitive root trait syndromes observed in plants growing in the wet sites (Fig. 3B). However, variation in root traits among sites was not significantly related with variation in bacterial community composition (P ≥ 0.609 for all seven traits; Table 2), while variation in bacterial communities did significantly correlate to two leaf physiological traits—mass based photosynthesis (P = 0.019) and water use efficiency (P = 0.018) (Table 2, Fig. 3C).

TABLE 2 Permutation tests performed on redundancy analyses predicting variation in bacterial and fungal communities as a function of environmental variables, root functional traits, and leaf functional traits^a

	Bacterial communities			Fungal communities
Traits	Variables	Variance	F value (P value)	Variables	Variance	F value (P value)
Environment	Canopy openness	0.017	1.86 (0.020)	Canopy openness	0.017	1.94 (0.009)
	Soil P	0.013	1.38 (0.106)	Soil P	0.016	1.79 (0.016)
	Soil pH	0.023	2.56 (0.001)	Soil pH	0.013	1.49 (0.051)
	Soil C	0.019	2.06 (0.011)	Soil C	0.016	1.86 (0.007)
	Soil N	0.011	1.15 (0.276)	Soil N	0.03	3.45 (0.001)
	Soil moisture	0.008	0.82 (0.664)	Soil moisture	0.016	1.86 (0.009)
	Residual	0.474		Residual	0.651
Root	Root diameter	0.015	1.55 (0.069)	Root diameter	0.015	1.67 (0.021)
	Specific root length	0.009	0.94 (0.506)	Specific root length	0.008	0.87 (0.662)
	Specific root area	0.007	0.76 (0.724)	Specific root area	0.009	0.94 (0.579)
	Specific root tip density	0.01	1.05 (0.389)	Specific root tip density	0.01	1.04 (0.381)
	Root length density	0.011	1.13 (0.299)	Root length density	0.013	1.39 (0.07)
	Root N	0.009	0.87 (0.601)	Root N	0.025	2.68 (0.001)
	Root C:N	0.008	0.79 (0.721)	Root C:N	0.014	1.55 (0.028)
	Residual	0.496		Residual	0.667
Leaf	Mass-based photosynthesis	0.019	2.03 (0.011)	Mass-based photosynthesis	2.25	0.02 (0.002)
	Water use efficiency	0.018	1.95 (0.014)	Water use efficiency	1.04	0.01 (0.364)
	Leaf density	0.014	1.46 (0.106)	Leaf density	1.59	0.02 (0.03)
	Leaf mass area	0.008	0.94 (0.508)	Leaf mass area	2.27	0.02 (0.003)
	Leaf N	0.008	0.81 (0.7)	Leaf N	0.9	0.01 (0.602)
	Leaf area	0.009	0.95 (0.468)	Leaf area	1.36	0.01 (0.104)
	Residual	0.488		Residual	0.67

^a Bacterial and fungal communities were assessed separately, with results based on 9,999 permutations. Significant predictors (where P ≤ 0.05) are highlighted in bold and results are presented graphically in Figure 3A to F.

View as image HTML

Download as PowerPoint

Variation in fungal communities across environmental gradients and with functional traits.

Fungal communities varied significantly as a function of all environmental parameters measured here including canopy openness, soil P, soil pH, soil C, soil N, and soil moisture (P ≤ 0.05 in all six cases) (Table 2). Compared with dry sites, wetter sites expressed higher values of soil P, C, N, and moisture, coupled with lower pH, all of which exerted a significant effect on fungal community composition (P ≤ 0.05 in all five cases) (Table 2, Fig. 3D). However, unlike bacterial communities, fungal community composition did vary with coffee root traits. Specifically, cool and dry conditions promoted the expression of high values of root diameter and root C:N, coupled with low root N concentrations; variation which was significantly related to fungal endophyte community composition (P ≤ 0.021 in all three cases) (Table 2, Fig. 3E). Fungal community composition also varied depending on two leaf morphological traits—leaf density (P = 0.03) and leaf mass per unit area (P = 0.003) (Table 2, Fig. 3F).

DISCUSSION

This study is the first report on root endophytes of coffee grown in Central America. Sequencing of all tree samples individually, rather than in pools, would have been preferable and even more informative; however, using a two pronged analytical approach, we were able to address two major kinds of questions. The TRFLP data allowed us to quantify microbiome variance within and between sites and to look at links to environmental factors and functional traits, but without identification of particular phylotypes responsible for differences. The sequence data on site/management pools allowed us to delineate organisms that appear to be core to Coffea arabica, and to note some that are putatively responsive to climatic factors.

We found 26 bacterial and 31 fungal genera that met our criteria as core—bacteria present at all four sites and average abundances above 0.5%, or fungi present at three or more sites and ranked in top 30 of average abundances for fungi (Table 1). While bacterial communities were consistent with the existence of a substantial core, the fungal community composition changes at the hot dry site are more suggestive of a gradient model (Hamady and Knight 2009). An even smaller number of taxa including 11 bacterial genera and 15 fungal genera comprised ≥40% of the total microbial relative abundance across a wide span of climate and environmental conditions (Fig. 2B and C). Relationships among these high-abundance microbial taxa varied widely, with certain genera being highly interconnected, and others being completely unrelated to other members of the coffee core microbiome (Fig. 2A). The most important result from this study is that coffee roots harbor both bacterial and fungal strains that are responsive to geography/climatic differences but not to the management regime (sun versus shade). This lends support to the concept that a core microbiome is conserved within coffee genotypes and only weakly influenced by contributions from surrounding soil and/or other plants.

There were some significant linkages between the microbiomes and some plant functional traits, particularly when looking at the fungi. It is most likely that both plant and microbes are responding to ecological selection imposed by the climatic differences, i.e., we are looking at covarying responses, as we have no evidence that the microflora are directly affecting plant functional traits. For instance, pH is a strong predictor of soil bacterial community compositions (Lauber et al. 2009), and soil bacterial communities are drivers of root bacterial communities (Bulgarelli et al. 2012). However, we cannot discount a direct effect of the microbiome on the ability of coffee to respond adaptively to changes in water and nutrient stress.

An inherent weakness of our data is that it is based on ribosomal operon-based sequence information only, i.e., phylogenetic data. Phylogeny is a very poor predictor of ecological function of both bacteria and fungi (Philippot et al. 2010). Ecophysiological functions can only be resolved by deeper, metagenomic studies or by reductive experiments on individual microbial isolates. The difficulty inherent in obtaining full genomic information on endophytic fungal or bacterial from DNA pools extracted from entire plants currently precludes metagenomic insights, and the equally inherent difficulty in culturing microbial associates without their hosts precludes reductive experiments. For now, we can only infer the importance of these microbes to coffee from their widespread abundance, and from their known associations with other plants.

Fungal members of the coffee root microbiome.

In the coffee trees we sampled, AMF were represented by Acaulospora, Archaespora, and Glomus (Table 1) with very minor amounts of Dentisculata and Funneliformis. Of these groups, Acaulospora and Glomus were on average the most abundant AMF, the other three were found at much lower abundances. Glomus and Acaulospora have also been found in coffee growing in more arid regions including Southwest Ethiopia (De Beenhouwer et al. 2015a) and Saudi Arabia (Mahdhi et al. 2017), as well as multiple other locations (Andrade et al. 2009). Here, our most striking observation was the complete absence of Acaulospora from the roots of coffee plants in both dry sites, and the absence of Glomus from the hot dry sun plantation (Table 1); patterns are consistent with AMF sensitivity to moisture stress and point to possible threats to plant−AMF associations under a shifting climate.

The most abundant fungal organisms isolated from the roots of these coffee trees were from genera traditionally perceived as pathogens or saprophytes. However, numerous examples of plant beneficial members of these genera can be found in the literature. Cladosporium made up about half of the fungal taxa found at the wetter sites. Cladosporium cladosporioides has been championed as a biocontrol agent against apple scab (Köhl et al. 2009) and other fungal pathogens (Wang et al. 2013). The majority of Cladosporium sequences found here were of this species. Hamayun et al. (2010) reported the production of gibberellic acid from a Cladosporium species isolated from cucumber root, and subsequent stimulation of cucumber shoot growth. Sette et al. (2006) isolated Cladosporium strains capable of antimicrobial production from coffee tissues. Similarly, Vega et al. (2006) note the isolation of numerous Penicillium strains from coffee tissues without proposing a function for them.

Numerous other non-mycorrhizal fungi with ambiguous ecological roles were also detected as part of the coffee microbiome. For example, both Colletotrichum and Paecilomyces were present in all sites and have also been observed in coffee growing in Mexico, Columbia, Hawaii, and Puerto Rico (Vega et al. 2010). Colletotrichum in particular is a large genus harboring well-known coffee pathogens, but also beneficial endophytes: Bongiorno et al. (2016) describe numerous Colletotrichum endophytes on coffee leaves in Brazil that possess antifungal activity, including activity against their congeneric pathogens. Hiruma et al. (2016) report root endophytic Colletotrichum tofieldiae strain that supplies phosphorus to Arabidopsis host in P-deficient soils.

The coffee endophytic microbiome contains several Ascomycota genera from orders known to harbor an entirely new and poorly understood group of fungal endophytes: the polyphyletic dark septate endophytes (DSE), which have been observed in a wide variety of terrestrial habitats. Neither their functions nor their taxonomy have yet to be fully elucidated (Mandyam and Jumpponen 2005). DSE in the coffee trees studied here include Codinaeopsis and Conoplea (Peziziomycotina) as well as Camarosporium and Stagonosporopsis (Pleosporales). While our understanding of the ecological role of many fungi is limited, many of these DSE can co-occur with AMF.

Two genera, namely Exidiopsis and Trechispora, were among the very few genera to exhibit pronounced preferences for hot and dry conditions (Table 1). Very little is known about either of these genera. The OTU we name as Exidiopsis had an ITS sequence that matched only poorly (90% nucleotide similarity) to any strain in the NCBI database. Similarly, Trechispora was also poorly matched to any known isolates or even nonculturable strains. However, given the higher relative abundances of these genera in hot and dry conditions, the ecological role of these groups could be expected to be important for understanding and predicting coffee−microbial associations under climate change.

The basidiomycete Mycena was found in roots at all sites in this study. This small, sometimes bioluminescent mushroom is thought of primarily as a colonist of decaying forest floor debris and senescing moss tissues, but there are numerous reports of Mycena acting as a mycorrhizal associate, providing carbon or other growth promoters to a wide range of different plant species and types (Grelet et al. 2017; Marquez et al. 2007; Zhang et al. 2012).

Also present in roots at all sites were Cryptococcus species. Some members of the genus are reported as beneficial endophytes of rice and sugarcane (Nutaratat et al. 2014), and in Brassica challenged with heavy metal pollution (Deng et al. 2012). Dokmaia, present only in roots at the wetter sites, is a novel genus only known from the description of saprophytes on the dead leaves of the endangered magnolia species Manglietia garrettii (Promputtha et al. 2003). There do not to seem to be any other reports of this genus in other systems.

Found abundantly at coffee growing at wet site Aquiares, and in lesser amounts at the other wet site, Trichoderma was represented largely by Trichoderma hamatum. A strain of this species was found to be highly effective against powdery mildew (Siddaiah et al. 2017). Another strain of the same species was found to protect Theobroma cacao seedlings during drought tests (Bae et al. 2009). High diversity of Trichoderma strains was also noted by Bongiorno et al. (2016) in organic coffee plantations.

Bacterial members of the coffee root microbiome.

The three most abundant bacterial genera that we found in coffee roots—Pantoea, Enterobacter, and Burkholderia—were nearly twice as abundant on average as any other bacterial genus (Table 1). Pantoea is a genus originally described to delineate a number of related agricultural pathogens found within the family Enterobacteriaceae, but later strains proved to be much more ecologically diverse (Walterson and Stavrinides 2015). Pantoea has been shown to be an ecologically important associate of many of the world’s most widespread crops, playing a role in growth promotion in winter wheat (Scholz-Seidel and Ruppel 1992), sugar cane, rice (Mano and Morisaki 2008), cut grass (Verma et al. 2018), some legumes, and rice (Megías et al. 2017), and as an N fixer in sweet potato (Asis and Adachi 2004).

Certain species of Enterobacter are critical to the growth of their host plants. Specifically, Enterobacter asburiae PDA134 (isolated from a date palm) carries genes for N fixation, auxin production, ACC deaminase, and siderophore production (Yaish 2016), while Enterobacter Sa187 isolated from a desert legume carries genes for siderophore production and iron update, and provides Arabidopsis some protection from salt stress (Lafi et al. 2017). Enterobacter cloacae MSR1 (isolated from Medicago plants) has been shown to enhance biomass gain in Pisum sativum (pea plants) (Khalifa et al. 2016).

Burkholderia, a large and very complex genus with large multireplicon genomes, contains species that exhibit close associations with multiple plant species (reviewed by Eberl and Vandamme 2016), and is now well known as a nodulating N₂-fixing plant in certain tropical plants (Walker et al. 2015). In addition to this, some Burkholderia strains form a unique kind of obligate symbiosis in the leaves of some African members of the families Rubiaceae and Primulaceae (Walker et al. 2015; Verstraete et al. 2017).

Two other genera well represented in coffee roots with likely important ecological roles in coffee are Erwinia and Pseudomonas. Two Erwinia and several Pseudomonas species isolated from coffee plants in Southwestern Ethiopia have been described as being highly efficient in P solubilization (Muleta et al. 2013), whereas multiple Erwinia strains have been found to be antagonistic toward plant-pathogenic fungi in crops (Misaghi and Donndelinger 1990; Li et al. 2010). Pseudomonas species are highly diverse, but numerous strains have been documented to have significant plant growth promotion abilities (Ma et al. 2017; Malfanova 2013; Mercado-Blanco et al. 2016; Padder et al. 2016; Pham et al. 2017).

In spite of it having the largest genome of any known prokaryote, nothing is known about the ecological role of Ktedonobacter (Chang et al. 2011; Yokota 2012). In this work, Ktedonobacter abundance increased, albeit slightly, under dry conditions. It also seems to act as a hub organism according to our network analysis and may play a role in mediating microbiome responses to the climate.

Conclusions.

Our data analyses here, when taken with studies from other parts of the world, provide support for the existence of a core microbiome of coffee. Our work shows that the bacterial members are much more consistent than fungal members of the core, but both vary across large environmental gradients and with intraspecific variation in plant functional traits. Surprisingly, the management regime (sun versus shade planting) had no significant effect on microbiome structure at any one site. Most genera found in the core microbiome contain well-known pathogens, but also several species and strains proven to have plant beneficial properties. If we are to understand the potential assistance that microbes can offer coffee and other crops in the face of our changing climate, in depth studies on the functional role of the key genera found here are suggested. It is notable that where beneficial associations involving these key coffee genera were noted in the literature, they were not restricted to coffee plants alone. In particular, we recommend study of the centrally connected Ktedonobacter and the poorly identified Exidiopsis, both of which were responsive to hot and dry conditions.

ACKNOWLEDGMENTS

We thank Bruno Rapidel, Olivier Roupsard, Karel Van den Meersche, Mirna Barrios, and Elias de Melo Virginio Filho for access to field sites and assistance with study design; and Shameek Das and Rosemary Saati for their invaluable lab assistance.