Phylogenetic Diversity of Arbuscular Mycorrhizal Fungal Communities Increases with Crop Age in Coffea arabica Plantations

Arbuscular mycorrhizal (AM) fungi are key soil microorganisms that establish a mutualistic symbiotic relationship with plants. The establishment of crops represents an environmental filter that usually reduces the diversity and variability of AM fungal communities, affecting the ecosystem stability and functionality. Despite several studies addressing these effects, the temporal development of these soil microbes since crop establishment has not been studied. We hypothesized that the negative effect of cropping practices in terms of reducing AM fungal richness, phylogenetic, and beta diversity will increase in time as far as the new dynamics progressively filter the AM fungal community composition. This research tested the impact of crop establishment and the role that time has in the progressive assembly of soil microbial communities. The AM fungal communities were characterized using terminal restriction fragment length polymorphism in coffee (Coffea arabica) plantations of different ages established in previous pristine tropical forest. We found that intraradical colonization and AM fungal phylogenetic diversity increased with plantation age. AM fungal richness was constant across time but a significant compositional turnover was detected. In relation to our initial hypothesis, these unexpected results face the current general view of the negative effects of crops on soil microbial diversity and highlight the need of studying temporal dynamics when assessing human impacts on soil biodiversity. Nevertheless, next steps would imply to put in context the found patterns by relativizing them to the original natural diversity inhabiting the studied areas.


Introduction
Soil microbial communities play a crucial role in terrestrial ecosystems due to their influence on different processes such as nutrient cycling or plant productivity (van der Heijden et al. 2008). This multifunctional aspect is related to the taxonomic and phylogenetic microbial diversity as far as functions are usually conserved in the phylogeny (Delgado-Baquerizo et al. 2016). Despite this acknowledge, the anthropogenic transformation of natural environments, specially by agriculture, often has a negative impact on soil microbial diversity (Jangid et al. 2008, but see García de León et al. 2018a, 2018b. Indeed, from an ecological perspective, crop establishment represents an environmental filter that causes the selection of particular microbial taxa as a consequence of the agricultural practices and the reduced availability of host plants (Heinen et al. 2020). In agreement, it has been observed that intensive agricultural practices, such as continuous crop monoculture, negatively affect soil microbial diversity (Mo et al. 2016), meanwhile alternative practices, such as, crop rotation, have been seen to increase soil microbial diversity (Venter et al. 2016).
Arbuscular mycorrhizal (AM) fungi are key components of soil microbiota. They establish a mutualistic symbiosis by which they colonize the roots of a major part of terrestrial plant species (van der Heijden et al. 2015). As mutualistic symbionts, they provide benefits to their host plants as increased nutrient uptake and protection against biotic and abiotic stresses (Powell and Rillig, 2018). Their presence and the composition of their communities affect plant diversity and productivity of ecosystems (van der Heijden et al. 1998). In turn, plant diversity and identity can also impact the diversity of associated AM fungal communities (Hausmann and Hawkes, 2009;Martínez-García et al. 2015). It is known that the effect of plant host identity on the AM fungal community is mediated by particular plant traits (López-García et al. 2017). Additionally, it has been shown that AM fungal community composition responds to the abiotic conditions of the environment (e.g., Chen et al. 2017) including disturbance (Trejo et al. 2016).
AM fungal communities are affected by land conversion to agriculture (e.g., Oehl et al. 2003;Verbruggen and Kiers, 2010). Their response to agriculture is variable and depends on several factors principally related to their functional spectra (Powell et al. 2009) or host plant traits. The output of the symbiosis between AM fungal taxa and particular host plants can also vary from positive to negative (Lehmann et al. 2012) and it has been observed that more aggressive agricultural practices tend to promote less mutualistic AM fungal traits, such as less quantities of extraradical mycelium or arbuscular colonization (Chagnon et al. 2022).
Despite the increasing evidence that points towards the general negative effects produced by intensive agriculture on AM fungal diversity and how they can be mitigated by adopting alternative practices, how the AM fungal communities respond to agriculture establishment in a temporal perspective (from establishment to subsequent crop seasons) has not been extensively addressed. This aspect is essential to gain knowledge on the response of microbial communities to initial disturbance that could determine the future composition and functionality in agricultural systems. Indeed, disturbances, such as the plowing and removal of previous vegetation implicit to land-use conversion into agriculture, have been shown to reduce overall AM fungal abundance (Oehl et al. 2003), species richness (Antunes et al. 2009) or caused drastic shifts in community composition (Violi et al. 2008;Schnoor et al. 2011). After a disturbance event, some studies have shown that the composition of AM fungal communities can change in time at individual plant root level (e.g., López-García et al. 2014). Due to the, a priori, more conservative management of perennial crops compared to the yearly restart of annual crops, potential recovery of AM fungal communities after crop establishment could be more evident in the former. However, temporal effects after agricultural conversion in perennial crops can be also confounded with the effect of the own ontology of host plants (Vukicevich et al. 2019).
The effect of crop establishment on soil microbial diversity, not only on AM fungi, has been typically approached by comparing values of alpha diversity (e.g., species richness, Newbold et al. 2015) across time or managements, but how beta diversity varies in time is not usually addressed despite that biotic homogenization caused by land-use change is of great concern (Gossner et al. 2016;Wang et al. 2022). Beta diversity is defined as the degree of compositional heterogeneity showed by a group of communities or samples in the same site. A decrease in values of beta diversity indicates a homogenization of community composition and hence a decrease in the local species pool with consequent implications for the stability of the ecosystem in the long term (Hooper et al. 2005).
To add light into how initial disturbance affects alpha and beta diversity of AM fungi during crop management, we focused in coffee plantations. Coffea arabica "coffee" is a perennial crop with a cycle of up to 20 years, of high economic value and with recognized global importance (Krishnan, 2017). It is a highly arbuscular mycotrophic plant (De Beenhouwer et al. 2015). Thus, coffee plants have a very reduced growth in the absence of mycorrhizal symbiosis and depend on this mutualism to ensure adequate nutrition, especially in eroded and low fertile soils as it is the case of tropical soils (Cogo et al. 2017). We hypothesize a decrease in alpha, beta, and phylogenetic diversity of AM fungi in time after initial disturbance caused by coffee establishment, suggesting that the environmental filtering caused by the new biotic and abiotic constraints would select an adapted group of AM fungal taxa (Verbruggen and Kiers 2010;Chagnon et al. 2022;Wang et al. 2022), i.e., those showing ruderal life-history traits (sensu Chagnon et al. 2013). Similarly, we expect a directional change in time of the community composition towards a common pool of AM fungal species adapted to anthropogenic systems (García de León et al. 2018b;Ohsowski et al. 2014) and particularly to coffee plant hosts. We formulated three specific objectives: (i) to reveal changes in the alpha and beta diversity of AM fungal communities in a temporal gradient; (ii) to study the turnover of the AM fungal community composition; and (iii) to evaluate potential functional consequences of crop establishment using the phylogeny of AM fungi as a proxy of their functionality (Powell et al. 2009;Chagnon et al. 2013).

Study Area
The study area was located at San Martín region in the northeastern zone of Peru (Fig. 1). This region is characterized by a humid and warm tropical and subtropical climate with a mean yearly rainfall of 1500 mm and temperature ranging between 23 and 27 °C. The altitude ranges between 139 and 3080 m.a.s.l. and includes pre-montane, lower montane, and montane humid and pluvial forests. Two seasons can be distinguished: a dry one, from June to September, and a rainy one, from October to May.

Sampling Design and Processing
The sampling was carried out during July and August 2018 during the dry season, in three sites of the San Martin region (Fig. 1, Table S1): Cashnahuasi (El Dorado province) (S 6°25′40.83″ W 76°52′28.41″-791 m.a.s.l.), El Mirador (Lamas province) (S 6°20′7.05″ W 76°33′24.33″-867 m.a.s.l.), and Shucshuyacu (Moyobamba province) (S 6°9′23.90″ W 76°53′48.01″-1027 m.a.s.l), which are characterized by their intense coffee production. In each site, three plantations of Coffea arabica L. differing in age were selected. Plantation age was retrieved from local land owners. Previous land-use varied across the selected plots: other crops such as mandarin, mango, maize-bean, or natural vegetation (Table S1). Four individual plants, separated at least by 5 m, were selected in each age class and site (n = 36). Plant position was georeferenced. Two subsamples per plant containing approximately 5 g of roots and 1 kg of soil were carefully collected between 5 and 25 cm depth from two equidistant points, each 0.5 m from the main stem of each plant. Roots were excavated to trace them back to the stem of the coffee plants to ensure they belong to the targeted plant. Subsamples were subsequently pooled, mixed and placed in polyethylene bags inside a cooler at 4 °C and then transported to Laboratorio de Biología y Genética Molecular in Universidad Nacional de San Martín (Tarapoto, Peru) and immediately processed. Samples were separated into roots and soil. Roots were washed, dried with paper, cut into 1-2 cm pieces, and homogenized in water. Aliquots of 500 mg were frozen at − 80 °C for molecular analysis. The remaining roots were kept in 70° ethanol for determination of AMF intraradical colonization. Soil samples were air dried at room temperature, sieved through a 2-mm mesh, and stored at 4 °C for a maximum of 2 weeks before further processing.

AM Fungal Intraradical Colonization
Approximately 2 g of roots for each plant were stained according to the methodology proposed by Vierheilig et al. (1998) using Parker Quink ink as colorant. Once stained, the roots were cut into 1-cm segments, mounted on slides, and examined in a compound microscope (20 ×) by intersection method (McGonigle et al. 1990).

AM Fungal Community Characterization
The molecular methods are explained in depth in Supplementary material Methods. Briefly, DNA was extracted from 100 mg fine roots using the cetyltrimethylammonium bromide (CTAB) protocol. The extracted DNA was subjected to terminal restriction fragment length polymorphism (T-RFLP) fingerprinting using AM fungal specific primers following the methodology described in López-García (2020). This allowed to obtain an operational taxonomic unit (OTU) × sample presence-absence matrix (TRFLP database approach), where the OTUs were delimited by the monophyletic clade approach (MCA, Lekberg et al. 2014), and a TRF peak profile abundance matrix based in the abundance of peaks produced by the two most variable restriction enzymes. In comparison with massive sequencing approaches, TRFLP has been seen to perform well for the non-highly diverse AM fungal group. Studies comparing TRFLP with massive sequencing have demonstrated that TRFLP is a good option to describe changes in AMF communities along environmental and temporal gradients (see, e.g., Varela-Cervero et al. 2016 andLópez-García et al. 2014, referring to massive sequencing results of Varela- Cervero et al. 2015 andLópez-García et al. 2017, respectively;and Martínez-García et al. 2015).
The OTU sequences found in the study are available in GenBank under the accession numbers MZ329082-MZ329146.

Statistical Analysis
The effects of age, site, and their interaction on soil properties, intraradical colonization, and OTU richness were tested by linear models. Homogeneity in the distribution of residuals was checked visually for model validation. Variables were square root, logarithm, arcsine, or BoxCox transformed according to the case. Spearman's correlations were performed to determine the relationships between intraradical colonization and soil properties.
The influence of age, site, and their interaction on the AM fungal community composition was analyzed by permutational multivariate ANOVA (PERMANOVA, using adonis function, vegan R package) where permutations were restricted to site (999 permutations). Jaccard distance for T-RFLP database matrix and Bray-Curtis distance for T-RFLP peak-profile matrix were used as measures of community dissimilarity. These distances were also used to carry out non-metric multidimensional scaling (NMDS) ordinations for visualization AM fungal community variation per age group across sites. Differences in multivariate dispersion, i.e., beta diversity, across plantation age per site were arranged using the linear models on the values obtained for each site using betadisper function (vegan R package).
The phylogenetic structure of the AM fungal community was addressed by calculating the standardized effect size of the mean pair-wise distance (SES-MPD) using the picante R package (Kembel et al. 2010). This index compares the mean phylogenetic distance between pairs of species in a sample with that obtained in a series of generated random communities. The index was calculated using a phylogenetic distance matrix between representative sequences of each OTU in the experiment (local pool) calculated using MEGA X (Kumar et al. 2018) (corrected using Kimura-2 parameters substitution model) plus the presence/absence matrix of OTUs. The independent swap algorithm was used to generate 999 randomized null communities to obtain the SES-MPD. SES-MPD values (z-scores) were used as dependent variable to test the effect of age, site, and their interaction via lineal modeling. Hence, the increasing or decreasing tendency of SES-MPD values in time would be interpreted as increasing or decreasing of phylogenetic diversity, respectively. To assess if particular site × age levels show significant phylogenetic clustering or segregation (negative or positive SES-MPD values, respectively, Webb et al. 2002), the mean value for each level combination was compared with the null distribution via t-test.
All statistical analyses were carried out using R 4.0.2 (R Development Core Team 2020).

Soil Properties
The effect of plantation age was significant for silt and cation Na + content (positive trends, Figure S1) (Table S2). Site, as random factor, seemed to indicate that some soil properties were different across studied sites (Table S3). For example, soil pH varied between sites: El Dorado had neutral soil (pH 6.5-7.2) whereas Lamas and Moyobamba showed strongly acidic soils (pH 3.5-5). All the soils were considered slightly saline (< 2 dS/m) and were characterized by having high percentages of organic matter content (> 4%) although some variation between the 3 sites were reported.

AM Fungal Community Composition
Nine OTUs were found in the study area using the monophyletic clade approach (Fig. 3), which corresponded to 75% of the theoretical total diversity ( Figure S2). Blast search against GenBank revealed that all OTUs belonged to phylum Glomeromycota and were grouped in two families: Glomeraceae (4) and Acaulosporaceae (5) (Table S4, Figure S3). Aca8 (identified as Acaulospora foveata) was found in 36.9% of sequenced clones followed by Aca6 (Acaulospora scrobiculata) which appear in 29.2% of the clones. Aca9 appeared in 15.4% and each of the rest of the OTUs (Glo1, Glo2, Glo3, Glo4, Aca5, Aca7) was found in less than 5% of the analyzed clones. OTU richness averaged 4.64 ± 0.36 and 7.56 ± 0.46 in database and peak profile approaches, respectively. It was not correlated with any soil property or affected by plantation age neither when using database nor peak profile approaches (Table S5).
The AM fungal OTU composition was affected by site and its interaction with plantation age in the database and the peak profile approaches (Table 1). The interaction between variables indicates that the effect of age is more evident in some sites than others. This result was graphically illustrated in the NMDS ordination of AM fungal communities (Fig. 4), where it seems that AM fungal  Table 1 Analysis of arbuscular mycorrhizal (AM) fungal community composition via permutational multivariate ANOVA. The effect of sampling site, plantation age, and their interaction was tested for AM fungal community obtained by terminal restriction fragment length polymorphism (TRFLP) in its database and peak profile versions. Jaccard and Bray-Curtis distances were used as measure of community dissimilarity for database and peak profile approaches, respectively. P values marked in bold indicate significant effect of predictors. R 2 indicates the proportion of explained variance in AM fungal community composition communities across time overlap more in El Dorado site than in Lamas and Moyobamba. Beta diversity was not affected by age, site, or their interaction (data not shown); thus, no convergence on community composition was detected among time and sites (Fig. 4).

Discussion
In one of the first studies approaching Peruvian Amazon region, we described AM fungal communities in terms of intraradical colonization, community composition, taxonomic, and phylogenetic diversity associated with coffee plantations in this region. As expected, due to the high mycotrophy of coffee, we found overall high levels of mycorrhizal intraradical colonization; however, this was highly depended on the nutritional level of the soil present in each sampled site. It is widely known that mycorrhizal colonization is affected by soil nutrient availability, especially phosphorus (Treseder 2004), that makes host plants reduce their investment in the AM symbiosis (Soudzilovskaia . The negative correlation found between soil nutrients and intraradical colonization ratified this assertion. Soil pH has a complex relationship and affects AM fungal community because it is directly related to the amount of available nutrients (Cruz-Paredes et al. 2021). Some previous studies have indicated a decrease in AM fungal colonization in acidic soils through the existence of positive correlations between soil pH and AM fungal root colonization (Palta et al. 2016;Liu et al. 2021); however, our results pointed out that AM fungal intraradical colonization exhibited a positive association with the exchangeable cations Al +3 + H + and a negative association with pH in the soil. These results could indicate that native AM fungal communities reported in our study are adapted to typical high acidic soils frequent in tropical and subtropical zones (Andrade et al. 2009). Similarly, soil texture also determines nutritional resources and water availability and consequently can influence the distribution and activity of soil microorganisms (Vieira et al. 2020). At this regard, we found a positive correlation between silt had and AM fungal colonization that could be explained by the active participation of silt in the formation of soil pores that improves water retention and air circulation generating a suitable substrate for AM fungal development (Banstola et al. 2021).
We also reported an increase in the intraradical colonization level over time. Hishi et al. (2017) found an increase in colonization with root age in Cryptomeria japonica plants and suggested that it could be due to the need for increasing the absorption efficiency that is lost during plant aging. In our case, despite that only fine and the least lignified roots were collected, the increase in the percentage of colonization was associated to an increase in the phylogenetic diversity of the AM fungal community in the older coffee plantations that was, in turn, driven by increased presence of members of Glomeraceae family. A high genetic variation of AM fungal, even within the same species, have been seen to affect important symbiotic parameters such as the colonization rate and growth of extraradical hyphae (Lee et al. 2013). This is especially true in case of an increase in members of Glomeraceae that are characterized by developing high degrees of intraradical colonization (Powell et al. 2009).
AM fungal total intraradical colonization also can be used as a proxy of AM fungal abundance in soil (Barceló et al. 2020). In this sense, an original degradation produced by crop establishment could, in one hand, reduce the colonization potential in the area, decreasing values of total Fig. 5 Standardized effect size of mean pairwise phylogenetic distance (SES-MPD) of arbuscular mycorrhizal fungal (AMF) communities in relation to plantation age mycorrhizal colonization, and, on the other hand, select for stress-tolerant AM fungi (Acaulosporaceae family, sensu Chagnon et al. 2013), with low presence into roots. In time, the gradual reestablishment of AM fungal populations would increase propagule potential and the presence of families more adapted to the agricultural lifestyle and with tendency to spread wider into roots, i.e., Glomeraceae members (Chagnon et al. 2013). It is important to note that the generalization of Glomeraceae as AM fungi with ruderal characteristic, although very widespread, is facing more recent findings showing a wide functional variation inside this family (e.g., Davison et al. 2021).
Our molecular analysis with TRFLP database showed nine OTUs belonging to Glomeraceae and Acaulosporaceae families. Despite the lack of previous reports in our study area, other studies using cloning-sequencing TRFLP approaches in crops have found similar number of AM fungal OTUs meanwhile revealed significant patterns across land-use managements (e.g., Wetzel et al., 2014). AM fungi belonging to Acaulosporaceae and Glomeraceae families have been previously seen in association with coffee (De Beenhouwer et al. 2014 and other crops in the Amazon area (De Assis et al. 2018;Stürmer and Siqueira, 2010). The broad distribution and high incidence of both genera in all of our study sites is faced with the general accepted view that most of ruderal AM fungal taxa (those with capacity to produce more spores in a shorter time, colonize from all inoculum types and a certain resistance to anthropogenic disturbance and stress) belong to Glomeraceae family but not to Acaulosporaceae (Verbruggen and Kiers, 2010;Chagnon et al. 2013; but see Davison et al. 2021). It is important to note the absence of other glomeromycotan families previously reported for coffee such as Gigasporaceae (Aldrich-Wolfe et al. 2020) whose presence is usually negatively affected by anthropogenic practices (De Souza et al. 2005;Chagnon et al. 2013). Nevertheless, we cannot discard that OTUs with reduced presence into roots, as Gigasporaceae (Varela-Cervero et al. 2015), could have not been revealed using this cloningsequencing and T-RFLP approach.
Our results showed that the composition of AM fungal communities was determined by the plantation age. By contrast, no differences, neither reduction nor increase, in alpha or beta diversity were observed in the temporal gradient. The change in the composition of AM fungal communities due to the influence of the plantation age have been reported in previous studies (Johnson et al. 2005;Kil et al.2014). It has been suggested that this variation could be caused by the different nutritional needs that host plants experience during their life stages and the differential functionality of AM fungal taxa that can be required by plants (Lu et al. 2019). However, the direction to which the composition of the AM fungal community changed in our study depended on the sampling site. In other words, both the youngest and the most mature AM fungal communities associated to coffee plantations were different between sites. This could be an effect of the variable soil properties that each sampling site presented which would determine the local pool of AM fungal taxa (Alguacil et al. 2016;Vieira et al. 2019). Indeed, the studied sites were highly variable in terms of, e.g., organic matter content, exchangeable cation content, or pH. Soil pH, for example, is a key factor that affects the composition of AM fungal communities at the local and regional level (Stürmer et al. 2018). It can have a direct influence by altering the physiological state of the AM fungi and indirectly by regulating the bioavailability of nutrients in the soil (Xu et al.2017). Hence, some AM fungal families have shown more or less pH tolerance, e.g., members of Acaulosporaceae can be found mainly in acidic soils whereas members of Glomeraceae prefer alkaline or neutral substrates . Other soil factors, namely, P, N, K, and organic matter content, have shown a direct influence on AM fungal community by conditioning the presence of particular AM fungal species (Carballar-Hernández et al. 2017;Cui et al. 2016;Jansa et al. 2014).
In addition to the effect on the composition of AM fungal community, plantation age had also influenced their phylogenetic diversity. Several studies have pointed out that phylogenetic diversity is related to the functional diversity of AM fungi and, hence, with their contribution to ecosystem services, such as the increase in the ability of survival of their host plants (Powell et al. 2009;Lee et al. 2013). Our results evidence a tendency to increase phylogenetic diversity with increasing plantation age what was contrary to our expectations: the selection of certain phylogenetic groups by coffee plants. Lazo (2012) reported that the coffee crop in Peru presents some limitations in terms of management. In our study area, crop management practices are usually limited at the time of establishment and decrease in frequency over time. This fact implies a reduction in disturbance in time from the initial coffee establishment and hence could revert the initial filtering towards particular adapted AM fungal taxa. Nevertheless, any of the found values of SES-MPD differed from null expectations, indicating that, despite the increasing tendency, the community assembly of AM fungi in our study was mostly determined by neutral, stochastic, processes (Webb et al. 2002).
Finally, the influence of site, plantation age, and their interaction on AM fungal community composition was found more clearly when using the TRFLP peak profile approach in contrast with the TRFLP database approach. TRFLP peak profile approach has the advantage of providing abundance data which allows improving the characterization of the community since it contributes to the quantitative component to the study (see, e.g., Mummey et al. 2005). This quantitative component in ecological studies is important because it is known that the abundance of some common species, and not the richness, lead the response of communities to their environment (Winfree et al. 2015).

Conclusions
As far as we know, this research is the first study that performs an exploration at the molecular level of arbuscular mycorrhizal fungi associated with a crop in Peruvian Amazon. This research revealed changes in arbuscular mycorrhizal fungal community composition after coffee establishment and, more interestingly, an increase in their phylogenetic diversity with time since establishment. These results directly incise in the knowledge about the sustainability of coffee cultivation due to the potential to foster the soil biodiversity. Further studies are needed for, for example, investigate the level of phylogenetic diversity in adjacent natural areas and compare with the coffee agroecosystem, or cover the annual natural variability by, e.g., including different seasons in the study designs.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Data Availability Raw data of this study will be deposited in Dryad repository upon acceptance.

Conflict of Interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.