The Coyhaique National Reserve (45°31′42.96”S, 72°1′51.95”W) is located in the Aysén Region of Southern Chile. In this reserve, deciduous forests with a predominance of Nothofagus pumilio constitute the typical vegetation. The native mature forests were affected by several wildfires in the region between 1930 and 1950 (Quintanilla Pérez, 2007; Fajardo and Gundale, 2015), after which natural regeneration resulted in native second-growth forests. Also, to protect the soil from erosion, plantations of various kinds of pines (mainly Pinus contorta) were part of afforestation programs organized by local authorities. However, the study area is not currently exposed to regular management practices since it is a protected area, and then the current forest community is the result of long-standing colonization.

Five ~200 m transects away ~300 m from each other (S1 to S5) were selected along a gradient in the forest cover composition present in the reserve trails (Supplementary Figure S1).

Forest cover at each site was estimated using the middle of the intercept point (Goodall, 1952) and the point-centered quarter (Mitchell, 2010) methods. The first one was used to calculate the relative coverage (RC) of each species in 2 m height, according to the formula RCi=niNT, where n_i is the frequency of the species i and N_T is the number of sampling points per site (Supplementary Table S1). The second method was used to calculate various indices from the field-collected data for each site (Supplementary Tables S2–S6). In this method, a reference point was selected every 20 m and divided into four quadrants. The closest tree was described at each of these quadrants by three parameters: species, distance to the reference point (d_i), and circumference at breast height (CBH). Then, three basic values were calculated: the radius of the tree trunk (r), according to the formula ri=CBHi2·π; the basal area covered by each tree with the formula Ai=π·ri2; and the total density (TD) of each species following TDi=1(diQ)2, where Q is the total number of sampled trees. After this, the relative density (RD) of each species was calculated with the formula RDi=niQ, where n_i corresponds to the total number of trees of each species. The absolute density (D) was calculated dividing the relative density by the total density according to Di=RDiTDi. Additionally, five other indices were calculated by the point-centered quarter method (Supplementary Table S7). The frequency (F) was calculated according to the formula Fi=jik, where j_i is the total number of points of each species, and k is the number of sampled points. Then, the relative frequency (RF) was calculated by dividing the F value of each species by the sum for all three species, according to RFi=Fi∑F. The coverage (C) of each species was calculated with the formula Ci=Ai·Dini, and the relative coverage (RC) by RCi=Ci∑C. Finally, the importance value (IV) was calculated according to the formula IVi=RFi+RCi+RDi, where the relative density (RD) of each species was calculated with the formula RDi=niQ, where n_i corresponds to the total number of trees of each species and Q is the total number of sampled trees. The importance value was relativized as Relative Importance (proportion of the total importance value) as it is easier to compare.

From the transect defined at each study site, five Peltigera frigida thalli and their associated substrates (Supplementary Table S8), alongside five soil samples from the surroundings (Supplementary Table S9), were collected. The lichen-associated samples were collected at least 1 m from the next closest to minimize the resampling of the same genetic individual. Besides, the soil samples were collected at least 1 m away from any lichen. The samples were placed in paper bags and transported in cooled containers. In the laboratory, the thalli were separated from the attached substrate with sterile brush and spatula. Thallus samples were stored in paper bags at room temperature, while the substrate and soil samples were sieved and stored in plastic tubes at −20°C.

From each soil sample, physicochemical parameters were determined as follows: pH was measured by potentiometry using a 1:2 ratio in 1 M KCl. Moisture content (MC; Steubing et al., 2002) and total organic matter (TOM; Sadzawka et al., 2006) were calculated gravimetrically before and after desiccation and calcination, respectively. Nitrate (N-NO₃⁻), ammonium (N-NH₄⁺), and phosphorus (P Bray) concentrations were measured spectrophotometrically from 1:10 soil extracts in deionized water, 1 M KCl, and Bray 1 extract, respectively (Bray and Kurtz, 1945; Steubing et al., 2002). Phosphorus sequential fractionation was determined by following Hedley’s protocol (Hedley et al., 1982) with the modification of do Nascimento et al. (2015). The soil was sequentially extracted with distilled water, 0.5 M NaHCO₃ (pH 8), 0.1 M NaOH, and mineral phosphorus (Pm) was calorimetrically quantified (Murphy and Riley, 1962). Subsequently, all samples were digested with H₂SO₄ and K₂S₂O₈ to determine organic phosphorus (Po) as the difference between total and mineral phosphorus (Po = Pt-Pm).

DNA was extracted from thallus, substrate, and soil samples using the DNeasy PowerSoil Kit (Qiagen), with 0.25 g per soil and substrate samples, as suggested by the manufacturer, and 0.15 g per thallus sample according to pre-experiments performed (data not shown). Modifications compared to the standard protocol included an overnight incubation at 4°C after adding the precipitation solution C2 (step 8) and an increase in time from 30 s to 1 min for the centrifugation steps after bead beating (step 6), ethanol washing (step 16), and elution (step 21).

From the DNA of lichen thalli, mycobionts were identified by analyzing the fungal 28S rDNA region amplified with primers LIC24R (Miadlikowska and Lutzoni, 2000) and LR7 (Vilgalys and Hester, 1990), and the ITS region amplified with primers ITS1F (Gardes and Bruns, 1993) and ITS4 (White et al., 1990). In addition, photobionts were identified using the cyanobacterial 16S rDNA region with the primers PCR1 and PCR18 (Wilmotte et al., 1993). All amplicons were sequenced with the forward primers using the Genetic Analyzer 3730XL (Applied Biosystems) using a sequencing service (Macrogen, Seoul, Korea).

Mycobiont identity was confirmed by a de novo phylogenetic analysis of the 28S and ITS sequences, carried out on the T-BAS web platform using the Peltigera reference tree (Carbone et al., 2019). In the case of the cyanobionts, the evolutionary history was inferred using the Neighbor-Joining method. For this, evolutionary distances were calculated using Kimura’s 2-parameter method, as suggested by MEGA7 (Kumar et al., 2016). The analysis involved 55 nucleotide sequences, including those from this work and reference sequences obtained from both cyanolichens and free-living cyanobacteria. All positions containing deletions were removed. Evolutionary analyses were performed in the software MEGA7 (Kumar et al., 2016). Finally, the phylogenetic tree was edited in iTOL (Letunic and Bork, 2019).

Using the DNA obtained from the thallus, substrate, and soil samples, quantification of the bacterial markers for phosphate solubilization, including the quinoprotein glucose dehydrogenase gene (gcd), alkaline phosphatase gene (phoD), acid phosphatase gene (phoN), phytase gene (appA), and phosphonoacetaldehyde hydrolase gene (phnX), was performed using quantitative PCR (qPCR) on an ABI 7300 Real-time PCR System (Kurth et al., 2020) and specific primers for each gene (Bergkemper et al., 2016a). These primers were reported as specific for bacteria, except phoD, which amplified genes related to some fungi but not Peltigera (Bergkemper et al., 2016a). Therefore, hereinafter we refer to phosphate solubilizing bacteria (PSB), although we cannot discard the presence of fungal genes related to P turnover. In addition, the 16S rRNA gene was also quantified as a marker for bacterial abundance (Kurth et al., 2020).

After each qPCR run, a melting curve analysis was performed to verify the amplicon specificity. The quantification of the target gene was conducted by using serial dilutions of plasmid-encoding the target genes (10⁷ to 10¹ gene copies μl⁻¹). Details on the used standards were published previously (Bergkemper et al., 2016b). All samples were diluted 1:4 to avoid inhibition during amplification caused by co-extracted humic substances based on a previous dilution test (data not shown). The efficiency of the qPCR was calculated as Efficiency (%) = [10(−1/slope) − 1]. Values were always above 85%, while the R² of the standard curve was above 0.99.

The abundance data of phosphate solubilization genes were log10 transformed and subjected to the Shapiro–Wilk test. The abundance of genes and physicochemical parameters were compared using one-way ANOVA and Tukey’s multiple comparison test. Non-metric multidimensional scaling (NMDS) analysis and pairwise Adonis test were performed based on Bray-Curtis distance, while redundancy analysis (RDA) was based on Euclidean distance and a stepwise selection of environmental variables using the ordiR2step function. All statistical tests and graphics were performed in R Studio (R version 3.5.3) using the vegan (Oksanen et al., 2017) and ggplot (Wickham, 2016) packages.

For this study, we selected five sampling sites representing different tree cover types. The first transect is located in a native N. pumilio mature forest (S1); the second one is spanning the transition between the native mature forest and a native second-growth forest (S2); the third one is located in the second-growth forest (S3); the fourth one is in the transition between the native second-growth forest and a second-growth forest of N. pumilio with the presence of pine specimens (S4); and finally, the last transect is the native second-growth forest with the presence of pine specimens (S5).

The first method used to estimate the forest cover, the middle of the intercept point method, calculates the relative coverage of the canopy of a tree species at a certain height (in this case, 2 m). The sum of this variable can reach values greater than 1 since the canopies of the trees overlap. According to this method, N. pumilio had the highest relative tree canopy coverage in all five sites. The sites S2 and S3, with second-growth forests, showed a less relative coverage of the tree canopy of N. dombeyi. Finally, the presence of P. contorta displaced the relative coverage of N. dombeyi in sites S4 and S5 (Table 1). The second method, the point-centered quarter method, calculates relative importance as the sum of the relative density, relative frequency, and relative volume for each species, being a way of measuring tree dominance. According to this method, N. pumilio is the most dominant tree species in all the forest sites, but it decreases in S4 with a concomitant increase of N. dombeyi, while P. contorta increased its dominance from S4 to S5 (Table 1). In summary, from S1 to S5, there is a decrease in the coverage of N. pumilio and an increase in that of P. contorta. N. dombeyi was also present but with no apparent coverage gradient among the sites (Table 1).

The physicochemical characteristics of the soils of the five sites are summarized in Table 2. Total organic matter, P Bray, soluble Pm (Pm-H₂O), labile Pm (Pm-HCO₃⁻), and labile Po (Po-HCO₃⁻) were significantly higher in S5 than in S1. Notably, soluble Pm was about 15-fold higher, while P-Bray was 3-fold higher in S1 than S5. Conversely, moisture content and N-NO₃⁻ tend to significantly decrease from S1 to S5. Finally, pH was higher in S2 and S4 than in S3 and S5, while N-NH₄ and Po-NaOH were similar in all sites.

Five specimens of P. frigida thalli were selected from each of the study sites. The identification of the mycobionts was confirmed by a phylogenetic analysis based on the 28S rRNA gene (Supplementary Figure S2) and the ITS hypervariable region (ITS-HR; Miadlikowska et al., 2003; Supplementary Figure S3). Further, all specimens selected were associated with the same photobiont haplotype namely Nostoc sp. C01 (Zúñiga et al., 2015) according to a 16S rRNA gene-based phylogenetic analysis (Supplementary Figure S4).

From the samples of the three different microenvironments (i.e., P. frigida thalli, substrates, and the surrounding soils) collected from the five sites (i.e., S1 to S5), we measured by qPCR the relative abundance of one gene that codes for an enzyme that triggers mineral phosphate solubilization (gcd), four genes that codes for enzymes involved in organic phosphate solubilization (phoD, phoN, appA, phnX), and the 16S rRNA gene as a bacterial marker (Figure 1). The quinoprotein glucose dehydrogenase gene (gcd; Figure 1A) was the most abundant marker in all soil samples, and unexpectedly its abundance was ~40,000 times higher in soils than in substrates and ~ 3,000 times higher in soils than in thalli in all the study sites. However, there were no significant differences in the abundance of this marker gene comparing the same microenvironment among the sampling sites. The abundance of the alkaline phosphatase marker gene (phoD; Figure 1B) also did not show differences when comparing the sampling sites, and although in the microenvironments of the extreme sites (i.e., S1 and S5) its abundance did not show significant differences, in the soil samples of the intermediate sites its abundance was significantly lower than in substrates and thalli. On the contrary, the abundance of the acid phosphatase marker gene (phoN; Figure 1C) was significantly less abundant in the extreme sites (i.e., S1 and S5) than in the other sites, independent of the microenvironment. The abundance of the phytase marker gene (appA; Figure 1D) and the phosphonoacetaldehyde hydrolase marker gene (phnX; Figure 1E) had the highest variation in soils throughout the sites, but appA showed the highest abundance in thalli and phnX in substrates. Finally, the 16S rRNA marker gene (Figure 1F) abundance almost did not change among sites, except for the soil samples from S5 that presented the highest abundance for this microenvironment. However, the abundance of the bacterial marker gene was significantly higher in thalli than in soils and substrates from the intermediate sites (i.e., S2 to S4).

NMDS multivariate analysis of the diversity of phosphate solubilization marker genes shows a clear separation of soil samples from those of lichen-associated microenvironments (i.e., substrates and thalli) in all the sampling sites (Figure 2A; Table 3). On the contrary, the separation of substrate and thallus samples, although statistically significative, is partially overlapped, being more evident in the intermediate than in the extreme sampling sites (i.e., S1 and S5). When comparing the diversity in the different sampling sites (Figure 2B; Table 4), there are two well-defined groups, the first one including S2, S3 and S4, and the second one including S1 and S5. Notoriously, the 95% ellipse encompassing the points of S1 is the largest one and encloses the 95% ellipse of S5, which is noticeably smaller.

To determine the effect of the edaphic and tree cover parameters on the abundance of phosphate solubilizing genes, we performed redundancy analyses (RDA) and ordistep for selecting the variables that significantly contribute to the model. The RDA analysis of soil samples, which explained 60.1% of the total variance, corroborates that sites S2, S3, and S4 grouped closely together and that S1 is the most diverse of all sites, in contrast to S5, which is the less diverse (Supplementary Figure S5A). The variables that significantly contributed to this model were the relative abundance of tree canopies of N. pumilio, N dombeyi, and P. contorta, pH, and the labile fraction of inorganic mineral phosphorus (Pm-HCO₃⁻; Table 5). On the other hand, the variance explained by the RDA analyses for substrates (Supplementary Figure S5B) and thalli (Supplementary Figure S5C) are not robust enough (36.7 and 47.5% of the total variance explained, respectively) to appropriately assert patterns of phosphate solubilizing bacteria diversity. In these cases, the only significant explanatory variables were the relative abundance of N. pumilio tree canopies and the labile fraction of mineral phosphorus (Pm-HCO₃⁻), respectively (Table 5).

Changes in the soil P content are influenced by multiple and diverse factors. One of them are wildfires, which in Patagonian forests has been reported to have long-term effects on soil properties, such as decreases in organic C and N and increases in pH and extractable P (Alauzis et al., 2004). Another factor is the influence of trees on soil chemical properties through their root exudates and the loss of their leaves in deciduous species (such as N. pumilio), contributing to changes in pH and the relative content and availability of macronutrients (Binkley and Valentine, 1991; Augusto et al., 2000; Frouz et al., 2009; Iovieno et al., 2010). Additionally, enzymatic activities of phosphorus solubilizing microorganisms could also affect the soil P content, which not only responds to the soil P availability but also to the balance with other nutrients, mainly nitrogen (Sorkau et al., 2018). In fact, Margalef et al. (2017) showed that the availability of nitrogen in soil is the main factor explaining phosphatase activity in temperate climates.

The quantification of the PSB marker genes in the soil samples evidenced two noticeable results. The first one was the high abundance of gcd, which codes an enzyme that allows the direct oxidation of glucose to gluconic acid, one of the main mechanisms of extracellular Pm solubilization studied in bacteria (Sharma et al., 2013; Sindhu et al., 2014; An and Moe, 2016). Therefore, the solubilization through gluconic acid of Pm from mineral complexes found in soils would be important in the sampling sites. There are other organic acids such as malic, lactic, citric and oxalic acids, which can also solubilize Pm (Sharma et al., 2013), but they have been less studied in this regard because the enzymes associated with their production play different intracellular roles and could be not relevant in the P cycle (Grafe et al., 2018). The second one was the pattern of abundance of the acid phosphatase marker gene (phoN), which was lower in the extreme sampling sites (i.e., S1 and S5) and more similar to each other than to the abundance in the intermediate sites. This could be related, in part, to the higher values of available nitrogen (ammonium plus nitrate) measured in the extreme sampling sites regarding the transition sampling sites (i.e., S2 and S4). In addition, the intermediate sites (S2-S3-S4) correspond to sectors of young trees and, therefore, they have a more active growth than those in the mature forest (S1) and then a greater demand for phosphorus that could explain the higher abundance of this marker gene. The measured pH in the soil samples (5.3 to 5.6), suggests that the activity of the enzyme coded by phoN could be more relevant than the activity of the alkaline phosphatase (encoded by phoD), as has been shown for other forest soils (Staddon et al., 1998; Herbien and Neal, 2008).

On the other hand, the growth dynamics of the tree species present in the sampling sites is different. P. contorta is a fast-growing exotic conifer, which would demand higher amounts of nutrients, such as P, compared to the native species N. pumilio (Fajardo and Gundale, 2015). Our results are in accordance with a comparative study about the P flux in coniferous and deciduous forests of the northern hemisphere, where the concentration of P in the forest floor under coniferous trees was about three times higher than in deciduous forests (Sohrt et al., 2017). Since the P content of coniferous foliage is larger than in deciduous one (Sohrt et al., 2017) and litterfall plays a critical role in nutrient cycling (Attiwill and Adams, 1993), the tree cover could affect the abundance of phosphate solubilizing microorganisms in soils. Indeed, in our study, the tree cover type strongly affected the abundance of PSB in the forest soil, being the relative coverage of all three tree canopies identified as a significant driver. The most noteworthy result regarding the effect of the tree cover influence on the PSB marker genes is the reduced diversity detected in the site with the highest presence of pines (S5) compared with the site of native mature forest (S1), although in the former the number of bacteria is higher.

Other significant explanatory variables were the pH, reported as one of the main factors controlling bacterial communities in soils (Fierer and Jackson, 2006), and labile Pm, which could be the primary source of P for PSB (Sims and Pierzynski, 2005). pH was higher in the transition zones (S2 and S4) than in the native second-growth forests independently of the presence of pines (S3 and S5), but we found a significant increasing gradient of total organic matter and various phosphorus fractions from the native mature forest (S1) to the native second-growth forest with the presence of pine specimens (S5). Although the phosphorus content values in our samples were in the range of soils in southern Chilean forests (Redel et al., 2008), differences were still observed between the sites.

It is important to mention that other factors, such as rhizospheric exudates and the association with mycorrhizae, could affect soil P dynamics (Chen et al., 2008). Interestingly, all the tree species of our study associate with ectomycorrhizal fungi (Bradbury et al., 1998; Godoy and Marín, 2019), which are highly efficient in P turnover by combining a high capacity to explore soils, the formation of polyphosphates, and the activity of acid phosphatase and phytase enzymes (Chen et al., 2008). However, they are affected by forest anthropogenic disturbances such as fires and afforestation (Marín et al., 2017), which could explain our observations that some bacterial markers related to the solubilization of organic phosphate (e.g., phoN) were lower in the soils of the sites where ectomycorrhizae would be contributing more importantly to the P turnover (i.e., S1 and S5).

Our results reinforce the concept that microorganisms are susceptible to changes in environmental conditions and management practices; thus, they are good indicators for assessing changes in soil quality or soil recovery from disturbances (Quintanilla Pérez, 2007), such as the afforestation with exotic trees after wildfires. A deep understanding of microbial communities in forest habitats with different forest cover is then essential for predicting the response of forest ecosystems to changes in environmental conditions (Lladó et al., 2017).

Thalli and substrates could be considered microenvironments with a high level of nutrient turnover (N, P or C), as biological soil crusts (Kurth et al., 2021) or the rhizosphere (Kuzyakov and Blagodatskaya, 2015). According to our results, the solubilization of Po from stable molecules such as phytate (through the enzyme coded by appA) or phosphonate (through the enzyme coded by phnX) seems to be more important in these microenvironments, respectively, than the solubilization of Pm, which would be more relevant in soils due to the higher abundance of the gcd marker gene. Phytate is accumulated by plants (Balaban et al., 2016), and to a lesser extent by microorganisms (Turner et al., 2021), as a P storage compound. The capacity to hydrolyse phytate has been reported as widespread among lichens, which potentially contributes to P capture from atmospheric deposits and plant leachates (Higgins and Crittenden, 2015). Phosphonates, on the other hand, can be found in membranes or excretions of microorganisms (Horsman and Zechel, 2017; Calcott et al., 2018) and groups of phosphonate biosynthetic genes have been reported in lichenized Nostoc (Gagunashvili and Andrésson, 2018). Therefore, considering the close contact between the thallus and the substrate in foliose lichens, as P. frigida, that the secondary metabolites produced by the lichen would diffuse to its substrate (Leiva et al., 2016), and that the main photobiont of Peltigera lichens were detected in both thalli and substrates (Zúñiga et al., 2017), the lichen-related microenvironments could be a rich medium in this type of compounds.

As an essential element, P plays a relevant role as a limiting macronutrient in the growth and activity of lichens. For example, it has been observed that a single immersion in phosphate solution can double the annual growth of the lichen Lobaria pulmonaria (McCune and Caldwell, 2009), as well as that fertilization with P, could increase nitrogen fixation in the cyanolichens Peltigera aphtosa and Peltigera polydactyla (Weiss et al., 2005). Considering the relevance of P in the development of lichens, different studies have been focused on the possible role that the lichen microbiota would be playing in functions related to phosphate solubilization processes. Indeed, several publications report potential PSB in lichen thalli (Bates et al., 2012; Grube et al., 2015; Sigurbjörnsdóttir et al., 2015). Thus, despite the lower abundance of organic phosphate solubilizers in comparison to mineral phosphate solubilizers, the former could play an essential role in recycling the old parts of the thallus (Grube et al., 2015), releasing phosphate that members of the lichen symbiosis could use. It is important to highlight that phosphatase and phytase activity have been reported in several species of lichens (Whiteford et al., 2010; Higgins and Crittenden, 2015), attributed to lichen enzymes being produced by the mycobiont. However, the primers chosen for this study were designed to amplify bacterial genes (Bergkemper et al., 2016a), which would indicate that the lichen microbiome is at least contributing to P turnover. However, since the lichens of this study are terricolous (i.e., their substrate is soil), the high abundance of the gcd gene, as the molecular marker of mineral phosphate solubilizers, is in accordance with analyses of metagenomes of soils where this is the most abundant gene of a set of 40 genes that encode enzymes related to the P cycle (Bergkemper et al., 2016b; Grafe et al., 2018).

Previous studies have found that thalli and substrates constitute microenvironments that differentially structure the associated microbial communities (Grube and Berg, 2009; Bates et al., 2011; Hodkinson et al., 2012; Ramírez-Fernández et al., 2014; Leiva et al., 2016). In particular, a recent study has shown that the bacterial community associated with P. frigida thalli is more than a mere extension of the microbiota of their substrates, with several bacterial being highly abundant in lichens but almost absent in substrates (Leiva et al., 2021). In our study, some markers of P turnover were more abundant in thalli, which could be related to the fact that part of bacteria related to this function would be recruited from the substrates and would later be enriched in the thalli. This proposition is related to Cardinale et al. (2012), who suggests that those microbial groups with specific functions within the microbiome would be acquired from the new environment colonized by lichen propagules because they are better adapted to local conditions.

In contrast to soils, PSB in substrates and thalli were only correlated with N. pumilio cover and labile Pm, respectively, and RDA analysis did not explain much of the variance. This supports our proposal that the microenvironments of lichens and their substrates act as an environmental buffer reducing the influence of forest cover composition on bacteria involved in P turnover, and that other factors driving the association with lichens could be of importance, like microbiome specificity.