Skip to main content

The influence of bioclimate on soil microbial communities of cork oak

Abstract

Background

Soil microbiomes are important to maintain soil processes in forests and confer protection to plants against abiotic and biotic stresses. These microbiomes can be affected by environmental changes. In this work, soil microbial communities from different cork oak Portuguese forests under different edaphoclimatic conditions were described by using a metabarcoding strategy targeting ITS2 and 16S barcodes.

Results

A total of 11,974 fungal and 12,010 bacterial amplicon sequence variants (ASVs) were obtained, revealing rich and diverse microbial communities associated with different cork oak forests. Bioclimate was described as the major factor influencing variability in these communities (or bioclimates/cork oak forest for fungal community), followed by boron and granulometry. Also, pH explained variation of fungal communities, while C:N ratio contributed to bacterial variation. Fungal and bacterial biomarker genera for specific bioclimates were described. Their co-occurrence network revealed the existence of a complex and delicate balance among microbial communities.

Conclusions

The findings revealed that bacterial communities are more likely to be affected by different edaphoclimatic conditions than fungal communities, also predicting a higher impact of climate change on bacterial communities. The integration of cork oak fungal and bacterial microbiota under different bioclimates could be further explored to provide information about useful interactions for increasing cork oak forest sustainability in a world subject to climate changes.

Peer Review reports

Background

Soil microbial communities are affected by environmental changes [1] and remain affected even after the end of extreme events, like drought [2]. Besides weakening plants, such extreme events have serious impacts on the structure of microbial communities [3], including the increase on the relative abundance of rare species [2] and emergence of pathogens able to infect plants [4]. On the other hand, stressed plants are able to recruit specific rhizosphere microbes, by changing the composition of their root exudates, thus modulating their adaptation and survival to such adverse conditions [5, 6]. Rhizosphere soil microbiome is composed by different microorganisms, including fungi and bacteria, which live in a complex environment [7] and respond differently to environmental changes [8,9,10,11]. Such microorganisms can establish symbiotic relationships with plants [12, 13] and interact with each other for influencing the biotic and abiotic environment [14]. For example, the interaction between mycorrhizal fungi and mycorrhiza helper bacteria (MHB) can increase root colonization and mycorrhiza formation, reducing plant stress [15, 16] and enhancing mineralization of organic phosphorus on soil [17]. Also, bacteria are important for the development of ectomycorrhizal fruit bodies [18]. The ecological services played by microbial soil communities, including plant growth promotion, stress tolerance, disease resistance, and their participation on biogeochemical cycles (mainly carbon, nitrogen and phosphorous) [19], turns the study of such soil communities an important issue for understanding climate change outcomes.

Climate models are showing that the Mediterranean region is a climate change hotspot [20]. The combined effects of decreasing precipitation and increasing temperature place Mediterranean forests as one of the most vulnerable ecosystems [21]. Among these, cork oak (Quercus suber L.) forests cover a large area of the western Mediterranean, having a high ecological and socio-economic importance [22]. The impact of climate changes on these forests have fostered cork oaks decline, which has been increasingly reported within the Mediterranean region [23]. Cork oak health is expected to be even more affected by drought events due to the emergence and spread of pathogens resulting from climate changes [23]. As generally recognized for woody plants [24], the microbial communities residing in cork oak forests could play a crucial role for preserving this ecosystem biodiversity and functionality. In soils of cork oak forest, the fungal communities have been better described than bacterial ones. Most studies relied on the identification of fungal communities using fruiting bodies surveys (e.g. [25]), culture-based methods (e.g. [26]) or root tips barcoding (e.g. [27, 28]). Fewer studies have used metabarcoding approaches for studying fungal communities associated with cork oak forests and the existent were focused on the fungal diversity found in different land-uses of these forests [29, 30]. For studying the impact of land management intensity on soil microbial diversity, a single study reported the simultaneous survey of bacterial and fungal microbiota by using a high throughput sequencing approach [31]. The integration of data obtained from fungal and bacterial microbiota is of utmost importance for understanding the functioning of any ecosystem and could provide important clues for the conservation of soil health, even when considering climate change scenarios [32]. Taking this into account, using a metabarcoding approach we aim to understand: (i) the divergence in microbial (bacteria and fungi) composition of cork oak forests located in regions with different bioclimates (combination of precipitation and temperature data), (ii) which edaphoclimatic factors have more influence in shaping these communities, and (iii) the correlation between fungal and bacterial biomarkers in different bioclimates. Final outcomes are expected to provide new information for a better management of cork oak forests under predicted climate changes.

Material and methods

Sampled forests

Soil samples were collected from eight cork oak Portuguese forests with different bioclimate classifications (Table S1; Fig. S1). Bioclimate classification was performed using Emberger climatogram and Emberger index Q = (2000 x P)/(M2 - m2), where P is the annual precipitation (mm), M is the maximum temperature in the warmest month (K) and m is the minimum temperature in the coldest month (K) [33,34,35]. Precipitation and temperature data spanning 10 years (2006–2016) previous to sampling were extracted from TerraClimate dataset [36], using QGIS [37] with a 0.05 buffer for each location. Sampled forests also presented different forest systems [sobreiral (about 400 trees/ha) or montado (60–100 trees/ha)], use of forest (wild forest or forest pasture) and tillage (tilled or non-tilled).

Cork oak soil samples

From each forest, cork oak trees [6 in Limãos (LI), Alcobaça (AL), Gavião (GV), Grândola (GR), Herdade da Contenda (HC-CT and HC-MA); and 5 in Parque da Peneda-Gerês (PG-ER and PG-RC)] were randomly selected and three equidistant soil core samples were collected from the middle of each tree canopy, at 2–10 cm in depth, in 2017 (Table S1). Soil samples were stored in sterile plastic bags at 4 °C, until processing. For physicochemical analyses, soil samples from each forest were thoroughly mixed and sent to a service provider (A2 - Análises Químicas, Portugal). Analyses comprised granulometry, pH (H2O and CaCl2), electric conductivity (μS/cm), organic matter, organic carbon, total nitrogen, carbon:nitrogen ratio and elements [phosphorus (P2O5), potassium (K2O), calcium (CaO), magnesium (MgO), sulfur, iron, manganese, boron, and sodium] (Table S2). The remaining soil was sieved (45-mesh sieve) to remove root tissues and other larger components of soil, thoroughly homogenized and three soil replicates per forest were created and stored at − 80 °C.

DNA samples preparation and Illumina sequencing of soil microbes

DNA extraction was performed for each soil replicate using DNeasy PowerSoil Kit (Qiagen, Germany), according to manufacturer recommendations. DNA amplification was assessed by PCR assays using a) ITS1-F (CTTGGTCATTTAGAGGAAGTAA, [38]) and ITS2 (GCTGCGTTCTTCATCGATGC, [38]) primers for the internal transcribed spacer 1 (ITS1) region of fungi, and b) 799F (AACMGGATTAGATACCCKG, [39]) and 1492R (GGTTACCTTGTTACGACTT, [39]) primers for the 16S V5-V9 bacterial region. The PCR reaction mixtures (25 μl) contained 1x Complete NH4 reaction buffer (Bioron GmbH, Germany), 200 μM of each dNTP (NZYTech, Portugal), 5 μM of each primer, 1.25 U DFS-Taq DNA Polymerase (Bioron GmbH, Germany) and 1 μl of DNA template (20 ng/μl). Amplifications were performed in a T100™ Thermal Cycler (Bio-Rad, USA), using the following protocol: initial denaturation for 4 min at 94 °C; 35 cycles of 30 s at 94 °C, 30 s at 52 °C (ITS1) or 54 °C (16S) and 30 s at 72 °C; final elongation at 72 °C for 10 min. PCR products were run on a 1% (w/v) agarose gel, stained with Green Safe Premium (NZYTech, Portugal). DNA samples that resulted in amplification were quantified using a fluorescent DNA quantification assay with dsDNA HS Assay Kit (ThermoFisher Scientific, USA) and detected with a Qubit 3.0 Fluorometer (ThermoFisher Scientific, USA).

DNA samples presenting a concentration higher than 20 ng/μl (three samples/forest) were sequenced using an Illumina MiSeq® sequencer with the V3 chemistry, through paired-end sequencing (2 × 300 bp) by a service provider (Genoinseq, Portugal). To determine fungal community, a first PCR reaction was performed to target hypervariable region of ITS2 using a pool of forward primers [ITS3NGS1_F (CATCGATGAAGAACGCAG), ITS3NGS2_F (CAACGATGAAGAACGCAG-3), ITS3NGS3_F (CACCGATGAAGAACGCAG), ITS3NGS4_F (CATCGATGAAGAACGTAG), ITS3NGS5_F (CATCGATGAAGAACGTGG-3), and ITS3NGS10_F (CATCGATGAAGAACGCTG)] and the reverse primer ITS4NGS001_R (TCCTSCGCTTATTGATATGC) [40]. Separately, to determine bacterial community a first PCR reaction was performed to target hypervariable V5-V7 region of 16S using the forward primer 799F-Y (AACMGGATTAGATACCCKG) and the reverse primer 1193R (ACGTCATCCCCACCTTCC) [39, 41]. A second PCR reaction, for each community, added indexes and sequencing adapters to both ends of the amplified target region, according to manufacturer’s recommendations [42].

Reads processing and data analysis

Raw reads were extracted from Illumina MiSeq® System in fastq format. Sequencing adapters and reads with less than 100 bases were removed with PRINSEQ version 0.20.4 [43]. Trimming, based on quality scores, was performed using default parameters in Sickle [44]. Before merging, correction of errors in reads was performed using Bayeshammer module from SPAdes package [45]. The merging of overlapping paired-end reads (Merged reads; Table S3) and further quality filtering was performed using Usearch version 8.0.1623 [46]. Fastq-mcf from ea-utils package [47] was used to filter merged reads based on sequence size. The software micca version 1.7.0 [48] was used to create a single FASTQ file (micca merge) and to discard sequences with an expected error rate greater than 1% (micca filter). This software (micca otu) was also used to cluster amplicon sequence variants (ASV), by using UNOISE3 protocol and chimeric sequences removal (Reads clustered into ASV; Table S3). Taxonomy was assigned to each ASV using a reference database [UNITE database version 8.2 [49] for fungi and qiime-compatible Silva release 132 [50] for bacteria] with micca classify. ASV assigned as unclassified were removed from bacterial dataset. In fungal dataset, those ASV assigned as unclassified, Viridiplantae, Algae or others not corresponding to fungal ASV were removed. Biom-format tables were created, and each dataset was subsampled using Qiime version 1.9.0 [51] to the number of sequences in the sample with the lowest number (Table S3; fungi to 26,787 that was found in GR1; bacteria to 25,087 sequences that was found in AL3).

Diversity and statistical analysis

All analyses were performed using the subsampled datasets in RStudio version 4.0.2 [52], except when stated otherwise. Fungal and bacterial richness (S) and diversity [Gini-Simpson’s index (1-D) and Shannon’s index (H′)] were determined for the different groups of samples using alpha() function of microbiome package [53]. While 1-D index measures the evenness of a community and represents the probability that two individuals randomly selected from a sample will belong to different species [54, 55], H′ index measures the diversity by taking into account the number of individuals as well as the number of taxa [54, 56]. These ecological parameters were compared between samples from different bioclimates, where hyper-humid comprised results from PG-ER and PG-RC forests, humid from LI and AL forests, sub-humid from GV and GR forests, and semi-arid from HC-CT and HC-MA forests. Statistical analysis were performed using stat_compare_means() from ggpubr package [57]. Rarefaction curves were computed to determine sampling effort using rarecurve() of vegan package. Fungal and bacterial abundance and richness were determined for the different bioclimates using microbiome [53] and phyloseq [58] packages.

Non-metric multidimensional scaling (NMDS) was performed using square root transformation of data for all samples to understand community distribution. NMDS was performed using vegdist() of vegan package [59] to calculate Bray-Curtis dissimilarity indices matrix, envfit() to understand goodness of fit (999 permutations) of environmental variables and metaMDS() of the same package to obtain the ordination graph. Goodness of fitness of the model was measured using Kruskal’s stress (values lower than 0.2 represent a good ordination; [60]. For evaluating differences on microbial communities from distinct bioclimates, ANOSIM was performed in Community Analysis Package 5.0 [61], using the Bray-Curtis dissimilarity matrix. This matrix was also used to perform multivariate statistical tests of significance (PERMANOVA), by using adonis() function of vegan package [59] with 999 permutations and homogeneity of group dispersions (PERMDISP). To determine which edaphoclimatic factors were more important in shaping fungal and bacterial communities a redundancy analysis (RDA) was performed using the package vegan [59]. Trend surface analysis was performed to transform latitude-longitude data into flat Cartesian coordinates using geoXY() of SoDA package [62] and polynomials of degree 3 were computed using poly() of STATS package [52]. Flat Cartesian coordinates (designated as cork oak forest) were included in RDA with other edaphoclimatic variables (bioclimate, forest system, tillage, and soil parameters). The function rda() was used to perform redundancy analysis, while anova.cca() was used to perform Monte Carlo permutation test (1000 permutations) and to evaluate the significance of global model. Most parsimonious model was found by forward selection of explanatory variables using ordistep() with 999 permutations. Multicollinearity of variables was evaluated using vif.cca() and excluded if variance inflation factors (vif) > 20. The variation explained by the model was tested using RsquareAdj() and Monte Carlo permutation test (1000 permutations) to evaluate significance. To understand the contribution of each variable for the variation of microbial community, variation partitioning was performed for the most parsimonious model using varpart() function. Statistical significance of variables and conditional variables was assessed using test_vp4() from comecol package [63].

The remaining analyses were performed using only ASVs classified up to genus level. The fungal and bacterial core community were related with bioclimate classification using heatmap.2() function of gplots package [64]. Linear discriminant analysis Effect Size (LEfSe) was used (LDA score > 4; p < 0.05) to discover bioclimate biomarkers [65], which consist in those microbial genera that could explain differences in bioclimate (the most significant environmental factor revealed by RDA). An adjacent matrix was created using Hmisc [66] and Matrix [67] packages, to perform a co-occurrence network of biomarker genera considering Pearson’s correlation coefficient higher than 0.75 (p < 0.05). Adjacent matrix was open in Gephi version 0.9.2 [68] to visualize co-occurrence network. The functional groups of biomarker fungal and bacterial genera were assessed using FunGuild [69] and FAPROTAX version 1.2.4 [70], respectively. For correlating biomarker genera and functional groups, correlation matrixes were computed using Pearson’s correlation coefficient with rcorr() function of Hmisc package [66]. Corrplot() function from corrplot package [71] was used to visualize correlation matrixes, by only presenting statistically significant correlations (p < 0.05).

Results and discussion

Fungal and bacterial dataset processing

Fungal and bacterial communities of soil from cork oak forests with different bioclimates were analysed by pair-end Illumina MiSeq platform, using ITS2 and 16S barcodes, respectively. Using ITS2 as a barcode, a total of 2,365,263 raw reads were obtained from 24 cork oak soil samples, ranging from 40,904 to 172,166 raw reads per sample (Table S3). Processed reads (1958803) contained around 3.2% of non-fungal ASVs, including 1.9% unclassified sequences, 1.3% Viridiplantae, 0.03% Alveolata and 0.008% Metazoa reads. Sequences corresponding to these ASVs were removed, resulting in 1,897,062 high quality fungal sequences, which were clustered into 11,997 ASVs (amplicon around 400 bp). To mitigate potential bias introduced by different sequencing depth in each sample, dataset subsampling was performed to the sample with least number of fungal sequences (GR1; 26,787 sequences). After subsampling, a total of 642,888 sequences were clustered into 11,974 fungal ASVs (Table S4A), which were assigned to 15 phyla, 37 classes, 91 orders, 170 families, 302 genera and 178 species. The same methodology was used for the bacterial dataset. From all cork oak soil samples (24), a total of 1,610,298 raw reads were obtained by using V5-V7 region of 16S region as a barcode, ranging from 50,673 to 91,233 raw reads per sample (Table S3). Processed reads (819012) contained 0.1% non-bacterial sequences, all of which corresponding to unclassified ASVs. After non-bacterial ASVs removal, a total of 818,092 high-quality bacterial sequences were obtained, which were clustered into 12,035 ASVs (amplicon around 400 bp). Bacterial dataset subsampling was performed to the sample with the least number of bacterial sequences (AL3; 25,087 sequences). After subsampling, a total of 602,088 sequences were clustered into 12,010 bacterial ASVs (Table S4B), which were assigned to 22 phyla, 53 classes, 107 orders, 142 families, 248 genera and 16 species. Although in this work the percentage of unclassified microorganisms was low (1.9 and 0.1%, when using ITS2 and 16S barcodes, respectively), there were many ASVs that did not match lower taxonomic levels, such as genus (Fungi 60%; Bacteria 46.9%) or species (Fungi 92%; Bacteria 99.7%). This result supports the assumption that much of the microbial diversity present in soils still remains to be discovered and studied [72, 73].

Cork oak forests microbial communities: diversity and composition

Microbial communities varied in richness and diversity in cork oak soils from different forests (Fig. S2) and bioclimates (Fig. 1). Rarefaction curves agreed with these results and revealed the good representation of microbial communities, especially when considering all soil samples from a certain bioclimate (Fig. S3). Although rarefaction curves for different individual forests did not reach a plateau (Fig. S3A, C), the sampling effort was good enough to represent microbial communities from different bioclimates (Fig. S3B, D). Taking this into account, and due to the meaningful ecological outcomes, results will be always discussed considering bioclimates. Soils taken from sub-humid bioclimates always presented the richest and more diverse microbial assemblages, both for bacterial and fungal communities (Fig. 1). Interestingly, semi-arid soils were particularly poor in fungal taxa and diversity (Fig. 1A) but were rich and diverse concerning bacterial communities (Fig. 1B). This result agrees with previous studies that described higher bacterial richness and abundance on soils subjected to drought [74, 75]. Such results have been explained by the theory of low pore connectivity, in which water-filled soil pores are poorly connected to each other in dried soils, affecting competitive interactions between bacterial species, which result in an increased bacterial diversity [74]. Previous studies also revealed the influence of climate on fungal communities [76]. In the present work, the most humid places were the ones that displayed higher fungal diversity and richness compared to semi-arid forests. This result could arise from the facilitation of fungal development under high relative humidity [77, 78]. Furthermore, this result is in line with the diversity found in ectomycorrhizae (ECM) root tips from cork oak stands [28], suggesting that was the higher diversity found in soils that contributed to the increased number of ECM root tips and not the promotion of mycorrhization process by itself.

Fig. 1
figure 1

Diversity of fungal (A) and bacterial (B) communities found in cork oak forest soils from each bioclimate. Each bioclimate is represented by two different forests [hyper-humid comprising PG-ER and PG-RC forests (n = 6); humid, LI and AL forests (n = 6); sub-humid, GV and GR forests (n = 6); semi-arid, HC-CT and HC-MA forests (n = 6)]. S represents microbial richness, 1-D Gini-Simpson’s index and H′ Shannon’s index. Different letters mean statistically significant differences, determined by ANOVA followed by Kruskal-Wallis test (p ≤ 0.05)

The dominant phyla in fungal community were Ascomycota and Basidiomycota, representing 41 and 22% of richness and 34 and 39% of abundance, respectively (Fig. S4). The prevalence of these phyla agrees with other studies in cork oak [29,30,31]. Agaricomycetes (21% richness; 37% abundance), Eurotiomycetes (10% each) and Sordariomycetes (9% each) were the richest and most abundant classes. Agaricomycetes was mostly represented by Agaricales (9% richness, 16% abundance) and Russulales (4% richness, 13% abundance). As expected from a Fagaceae forest [79], ectomycorrhizal species dominated in diversity (18% of richness), but particularly in abundance (37%), mainly regarding the fungal families Russulaceae (13% of total abundance), Inocybaceae (5%), Cortinariaceae (3%) and Thelephoraceae (3%). These results are in line with other fungal metabarcoding studies on Fagaceae forests that also described an enriched ECM Basidiomycota community [80, 81]. Besides, these ECM families were considered as ubiquitous in Mediterranean cork oak forests [82]. When comparing the composition of soil forest fungal assemblages among different bioclimates, a similar richness profile was found (Fig. S5A), but differences in fungal abundance were detected in different bioclimates (Fig. 2A). More humid forests were particularly enriched in Russulaceae [ECM, 16.2% in most humid (hyper-humid and humid) forests against 9.2% in driest (sub-humid and semi-arid) ones] and Mortierellaceae (non-ECM, 9.6% against 3.6%), but also displayed high levels of Amanitaceae (ECM, 4.5% against 0.6%) and Cantharellaceae (ECM, 1.5% against 0.2%). The opposite pattern was observed for Inocybaceae (ECM, 2.7% against 7.6%), Cortinariaceae (ECM, 1.6% against 4.8%), and Hydnangiaceae (ECM, 0.4% against 3%) families. The same family distribution trend was observed for ECM root tips [28], suggesting once more the relation of ECM fungal abundance and mycorrhization of cork oaks. Other non-ECM families also displayed higher abundance in less humid bioclimates, like Piskurozymaceae family (0.2% in most humid forests against 2% in driest ones). But the most noteworthy difference of non-ECM fungi among different bioclimates was the higher occurrence of the abundant Mortierellaceae family in humid and sub-humid forests. Besides being humid, these forests are characterized by lower annual temperatures throughout the year, which could favour the development of these fungi, as warmed soils were recently reported to reduce the abundance of Mortierellaceae members [83]. As the most humid forests are mostly wild forests with high density of trees, the prevalence of Mortierellaceae members in these forests could be also related with their ecological function as saprophytic fungi, decomposing forest litter through chitin and lignin degradation [84]. However, bioclimate seems to be more important for Mortierellaceae abundance (PERMANOVA: R2 = 0.66, p = 0.001) than forest use (PERMANOVA: R2 = 0.25, p = 0.001).

Fig. 2
figure 2

Mean relative abundance of soil fungal (A) and bacterial (B) communities in forests from each bioclimate. Results display communities at family level, being indicated the corresponding classes and phyla

Regarding bacterial community, Proteobacteria, Actinobacteria and Acidobacteria were the most prevalent phyla with 35, 29 and 19% richness and 37, 32 and 18% abundance, respectively (Fig. S6). Similar profiles have been detected in other forest soils (e.g., [85]) and cork oak forests ([75]). The richest and more abundant classes were Alphaproteobacteria (19% richness; 22% abundance), Actinobacteria (17%; 21%), Acidobacteriia (16%; 17%) and Gammaproteobacteria (12%; 13%). Of relevance was the dominance of Xanthobacteraceae (Alphaproteobacteria, 10% of total abundance), followed by Acidothermaceae (Actinobacteria, 6%) and Burkholderiaceae (Gammaproteobacteria, 5%) families, all including well-known members able to interact with plants. For example, Xanthobacteraceae comprise species with plant growth promotion abilities [86] and nitrogen fixation in legumes [87]. The highest prevalence of Acidothermus sp. (Acidothermaceae) and Burkholderia sp. (Burkholderiaceae) have been described in different forests, including cork oak forests [75, 88]. As detected for fungal communities, a similar profile of richness was found in cork oak forests from different bioclimates (Fig. S5B). Regarding taxa abundance, differences between bioclimates were perceived, but not as clearly as in fungal community (Fig. 2B). Besides presenting a general high abundance, Acidothermaceae and Burkholderiaceae, as well as Mycobacteriaceae were more abundant in the most humid forests (4.2, 6.6 and 4.8%, respectively) than in the drier ones (2.6, 4.2 and 2.6%, respectively). The opposite pattern was shown for Solirubrobacteraceae (2.9% against 4.8%) and Nitrosomonadaceae (1.8% against 3.1%). The most arid bioclimates also seem to favour bacterial proliferation of Geodermatophilaceae (0.4% in most humid against 3% in the driest forests). Taken together, these results suggest a bioclimate impact on the abundance of some microbial taxonomic groups.

Environmental factors influencing microbial community structure

The influence of different edaphoclimatic factors (Table S1 and S2) in shaping fungal and bacterial communities was evaluated by fitting the environmental factors onto a NMDS ordination plot (Fig. 3). NMDS analyses shown a good representation of the distribution of microbial communities (Kruskal’s stress < 0.2), revealing that they grouped differently according to the sampled forest, location and bioclimate. The microbial communities from semi-arid regions were the most distant from the other communities in moister regions. These results were corroborated by ANOSIM analyses that revealed highly dissimilar communities in different bioclimate groups (fungi: R = 0.954, p = 0.001; bacteria: R = 0.935, p = 0.001). These results were reinforced by homogeneity of group dispersions (fungi: F-value = 44.911; p = 4.52e− 09; bacteria: F-value = 12.281; p = 8.87e− 05). Less dissimilarities were found among humid and sub-humid forests (fungi: R = 0.783, p = 0.001; bacteria: R = 0.780, p = 0.001), as well as between sub-humid and semi-arid bioclimates (bacteria: R = 0.665, p = 0.001). Communities between all the other bioclimates were almost completely dissimilar (R ≥ 0.941, p = 0.001 for both microbial communities). These results suggest that particular drivers could have specifically shaped the microbial communities in arid and humid environments.

Fig. 3
figure 3

Evaluation of the main environmental factors affecting fungal (A) and bacterial (B) diversity. Non-metric multidimensional scaling (NMDS) analyses were performed using Bray-Curtis dissimilarity coefficients. Each symbol represents a different soil sample (three samples from each forest), with colours representing different forest bioclimates. Environmental variables were fitted onto the NMDS ordination. The arrows directions indicate positive correlations between continuous environmental factors and microbial assemblages. Tmax represents the maximum temperature in the warmest month, Tmin is the minimum temperature in the coldest month), Organic C is the organic carbon, Total N is the total nitrogen and C:N ratio is the carbon:nitrogen ratio. Only those factors that were significantly correlated with NMDS ordination axes (p < 0.05) are represented

According to the environmental factors fitted to NMDS ordination plot, most of the tested edaphoclimatic parameters influenced the microbial communities, with exception of ‘Tmin’, ‘forest system’, ‘tillage’, and ‘total N, P, S, Fe’ (Fig. 3). For both microbial communities, the highest R2 values were found for ‘cork oak forest’, ‘forest location’, ‘bioclimate’, and ‘precipitation’ factors (all above R2 > 0.8), but ‘Tmax’, ‘boron’ and ‘C:N ratio’ (in both communities) and ‘manganese’ (in fungal communities) were also particularly relevant (R2 > 0.7). The importance of bioclimate in shaping fungal and bacterial communities is in line with previous studies that revealed the influence of climate factors in bacteria [75] and ectomycorrhizal tips distribution in cork oak soils [28]. In the present work, the redundancy analysis (RDA) also revealed ‘bioclimate’ as the factor that better explained the differences among microbial communities (Fig. 4). A higher contribution of ‘bioclimate’ was found when considering the variation of bacterial communities (R2 = 0.903, p = 0.001, NMDS; 28.1% of variation, p = 0.001, RDA) in relation to variation of fungal communities (R2 = 0.889, p = 0.001, NMDS; 9.8% of variation, p = 0.001, RDA). However, ‘bioclimate’ and ‘cork oak forest’ were collinear variables (variables are correlated) for fungal community and these results should be taken with caution. Detected differences could be related with the different role of drier environments in shaping bacterial and fungal communities. While driest samples presented a significant reduction on fungal richness and diversity, bacterial communities increased both ecological parameters in semi-arid forests (Fig. 1). These results agree with other studies that revealed a much stronger impact of drought on bacterial than on fungal networks [89].

Fig. 4
figure 4

Contribution of edaphoclimatic factors in shaping fungal (A) and bacterial (B) communities. The most parsimonious model resulting from redundancy analysis is shown. The best explanatory variables were identified and their combined contribution for microbial shaping is displayed. Statistically significance values based on 999 permutations are represented as *** p < 0.001, ** p < 0.01 and *p < 0.05

Besides ‘bioclimate’, both microbial communities were strongly affected by ‘cork oak forest’ and distance among cork oak forests (‘forest location’), both at R2 > 0.947, p = 0.001 (NMDS). But, neither of these factors revealed to significantly contribute to the shaping of bacterial community (Fig. 4). For shaping fungal communities besides collinear variables ‘bioclimate’ or ‘cork oak forest’, RDA displayed soil ‘pH’, ‘boron’ and ‘granulometry’ as playing a significant role in fungal community shaping (7.1%, p = 0.001, 6.7 and 3.8%, p = 0.001, respectively). All the other environmental factors were found to not contribute significantly for shaping fungal communities. A different picture was revealed by the factors that mostly shape bacterial communities (Fig. 4). The ‘C:N ratio’ and ‘boron’ factors played a major role in shaping bacterial communities (21.9%, p = 0.001 and 18.6%, p = 0.001, respectively), followed by soil granulometry (9.4%, p = 0.001). Altogether, ‘bioclimate’ and these three factors explained more than 60% of bacterial variation. This contrasts with the significant factors that shape fungal communities (‘bioclimate/cork oak forest’, ‘pH’, ‘boron’ and ‘granulometry’) that together only explained 27.4% of fungal variation. These results suggest that bacterial and fungal communities are shaped by different edaphoclimatic features at distinct levels, as previously reported by Shen et al. [90]. These authors suggested that pH and temperature were the best predictors for bacterial and fungal composition, respectively. In the present work, communities were not significantly affected by temperature alone, but bioclimate was distinguished as the main driver of bacterial communities. Interestingly, soil pH only explained variations detected on fungal community and was not significant for bacterial composition. This contrasts with multiple studies that report pH as the most important driver for bacterial diversity (e.g. [90, 91]). Detected differences could be related with the influence of other soil factors for shaping bacterial community. On the other hand, pH was found to be a significant factor influencing arbuscular communities [92] and was considered as the second most important driver of fungal communities on cork oak forests in Morocco, as reported by Maghnia et al. [29]. In the present work, pH was also the second most relevant driver for structuring studied fungal communities. Other edaphic features were also relevant for driving microbial communities, like it was the case of granulometry (for both communities). Accordingly, distinct soil granulometry contributed to different fungal communities in cork oak forests [26] and impacted bacterial diversity [93]. In addition, the structure of microbial communities was conditioned by ‘boron’ (both communities) and ‘C:N ratio’ (bacterial community). Boron is considered as a micronutrient essential for plant growth, but in high doses is known to influence microbial community composition by inhibiting soil microbial respiration [94]. However, fungi have been reported as more sensitive to boron in relation to bacteria [95], which contrasts with the finding of the strongest influence of boron levels on driving bacterial assembling when compared to its effect on fungal community shaping. Also, C:N ratio (representing the soil nutrient availability) modulates the bacterial community, as detected in different soil systems (e.g. [96, 97]). Although in the present work different land management practices (forest use and tillage) revealed none or little direct impact on soil microbial communities, they may have an indirect influence due to their influence on different soil properties, such as pH, carbon, nitrogen and phosphorus [98, 99]. The prominent influence of ‘bioclimate’ for the structure of microbial communities dictates the urgent need to understand the impact of climate change on microbial communities and final outcomes for cork oak forests sustainability.

Correlation of microbial communities in different bioclimates

The microbial taxa present in all soil samples included 102 fungal genera and 188 bacterial genera (core communities), which presented different abundances in distinct forest bioclimates (Fig. S7). Among fungi, many core ECM genera (Russula, Inocybe, Laccaria, Cenococcum, Amanita, Lactarius, Tuber, Sebacina, Tomentella and Hebeloma) have been previously described as core genera of Mediterranean cork oak forests [82]. The abundance of genera from the core fungal community in hyper-humid and humid forests clustered together, differently from the cluster formed by genera abundance in sub-humid and semi-arid forests. In contrast, the abundance of core bacterial community found in sub-humid and semi-arid bioclimates were the most similar, different from those present in humid forests, and even more dissimilar from the ones present in hyper-humid forests. The most abundant core genera were Russula (Russulaceae, 22% of core reads), Inocybe (Inocybaceae, 10%), Mortierella (Mortierellaceae, 8%), Penicillium (Aspergillaceae, 7%) and Amanita (Amanitaceae, 5%) (Fig. 5A). Following the general trend (Fig. 2A), Russula, Mortierella and Amanita decreased their presence from most humid to most arid bioclimates, which was particularly evident in Russula genus. In contrast, Inocybe and Penicillium were more present in semi-arid and sub-humid forests, respectively. Regarding the core bacterial community, the most abundant core genera were Acidothermus (Acidothermaceae, 11% of core reads), Bradyrhizobium (Xanthobacteraceae, 10%), Mycobacterium (Mycobacteriaceae, 7%) and Burkholderia-Caballeronia-Paraburkholderia (Burkholderiaceae, 5%) (Fig. 5B). These genera were mostly correlated with the most humid bioclimates. The families of all these bacterial genera were already reported as core families in cork oak forests [75].

Fig. 5
figure 5

Heatmap depicting the most abundant core fungal (A), and bacterial (B) genera present in each bioclimate. Core communities were assessed considering those ASVs present in all bioclimates and classified up to genus level. Color represents differences in abundance, where dark red represents high abundance and dark blue low abundance. This figure is a detail of the heatmap depicting all core genera represented Fig. S7

For identifying biomarker taxa for each bioclimate, a linear discriminant analysis was performed with all dataset. Different fungal and bacterial genera revealed to be biomarkers for the analysed bioclimates at a LDA score higher than 4 (Fig. S8). A total of 49 fungal (10 hyper-humid; 18 humid; 11 sub-humid; 10 semi-arid; Fig. S8A) and 69 bacterial genera (10 hyper-humid; 23 humid; 10 sub-humid; 26 semi-arid: Fig. S8B) were defined as biomarkers with statistical significance (p < 0.05). Most of them (34 fungal and 66 bacterial taxa) were also identified as core genera, reinforcing the previous suggestion that specific bioclimates may favour the colonization by certain genera, even when they are generally widespread in all cork oak soils. As expected, correlations between fungal and bacterial biomarkers from a specific bioclimate were positively correlated, mainly concerning those from the most extreme (hyper-humid and semi-arid) bioclimates (Fig. 6; Fig. S9). Also, fungal biomarkers from the hyper-humid forests revealed the most negative correlations with bacterial biomarkers from other forests, particularly from semi-arid forests. Although not so evident the same trend was observed for bacterial in relation to fungal biomarkers. This result suggests that genera found in most humid (hyper-humid and humid) bioclimates are not prone to develop in the most arid (sub-humid and semi-arid) ones and vice versa. Interestingly, certain biomarker fungi revealed a complementary correlation with bacterial biomarkers from all the other bioclimate. For example, Clavulina (humid biomarker) was positively correlated with all biomarker bacteria from all bioclimates, except from those of hyper-humid. On the other hand, Lactifluus (humid biomarker) displayed positive correlation with most humid bacterial biomarkers and negative from those of drier forests. The same was found for certain biomarker bacteria (e.g., Burkholderia-Caballeronia-Paraburkholderia vs. Sphingomonas). These results suggest that bioclimate modulates a complex co-occurrence network within both communities, determining microbial interactions with possible outcomes to forest sustainability. The recognition of physical and molecular fungal-bacterial interactions for plant health and for ecosystem functioning has been reported [13, 14]. For example, fungal communities were recently described to have an effect on bacterial composition in deadwood [100].

Fig. 6
figure 6

Co-occurrence network of fungal and bacterial biomarkers of different bioclimates. Pearson’s correlation coefficient was calculated for the top 50 genus displayed in LEfSe (R2 = 0.75, p < 0.05). Fungal nodes are represented by dashed circles and bacterial nodes by full circles. Size of the node corresponds to node degree (bigger nodes represent higher number of connections). Nodes are colored to represent biomarkers from different bioclimates (Hyper-humid in blue, Humid in yellow, Sub-humid in orange and Semi-arid in brown). Weighted edges in blue represent positive correlations and red represents negative correlations. Details of Pearson correlation matrix can be found in Fig. S9

The ecological functions of identified fungal and bacterial biomarkers in the ecosystem were predicted using FunGuild and FAPROTAX, respectively, and the relative abundance in each bioclimate is displayed (Fig. 7). The main ecological functions of biomarker fungi were as ectomycorrhizal (ECM) and saprophytic (SPT) fungi, being ECM prevalent in the most extreme environments (hyper-humid and semi-arid bioclimate, Fig. 7A). Indeed, ECM distribution have been related with climatic conditions (e.g., [101]), including in cork oak forests [27, 28]. The identified bacterial biomarkers were mainly related with chemoheterotrophy (33.30%) and aerobic chemoheterotrophy (32.97%). Bacteria with these ecological functions presented significantly higher abundance in hyper-humid than in semi-arid bioclimates (Fig. 7B; Table S5). Also, bacteria that perform cellulolysis represented 6.63% of total abundance and presented significantly higher abundance in hyper-humid bioclimate compared with all other bioclimates. These results suggest that different functional bacterial groups are thriving in different bioclimates, being the hyper-humid bioclimates more related with plant matter decomposition. When correlating different ecological functions within microbial communities (Fig. 7C), a significant negative correlation was observed among ECM and SPT fungi. Besides, these fungal guilds correlated differently with all the other bacterial functions. The different lifestyle of ECM and SPT and distinct nutrient acquisition strategies are able to create soil niches with discrete nutrient provisions that allow the propagation of bacteria with distinctive metabolic features [102]. Other ECM have been highly correlated with well-known mycorrhizae helper bacteria, and also with some not so well described bacteria, which reinforces the importance of fungi-bacteria interactions for forest sustainability [16, 103]. These results suggest once more the complex and delicate balance that occurs among microbial communities that is strongly dependent on bioclimate.

Fig. 7
figure 7

Relative abundance of functional groups of fungal (A) and bacterial (B) biomarkers in the different bioclimates and Pearson’s correlation matrix between functional groups (C). Only fungal and bacterial biomarkers obtained by LEfSe were used. Only more abundant bacterial functional groups are displayed in B. Full list of bacterial functional groups can be found in Table S5. ECM refers to ectomycorrhizae, SPT to saprotrophs and ERM to ericoid mycorrhiza. Blue and red circles represent statistically significant positive and negative correlations at p < 0.05, respectively. Color intensity is proportional to the correlation coefficients and circle size to statistical significance

Conclusions

Soil microbiomes are essential to maintain forest health and sustainability. Their role will be even more important under the pressure of ongoing climate changes. Among others (pH, boron levels, granulometry for fungi; C:N ratio, boron levels, granulometry for bacteria), bioclimate was the factor that contributed the most to shaping bacterial communities, while bioclimate/cork oak forest influences the structure of fungal communities in cork oak forests. These reinforces the concerns about these forests’ sustainability under a climate changing scenario. Forest soils with less water availability can develop richer and more abundant bacterial communities in comparison to fungal communities. However, differences on specific taxa abundance among bioclimates are more notable in fungi than in bacteria, raising the question about the effects of bioclimate in both microbial communities. The higher contribution of bioclimate for shaping bacterial communities predicts a higher impact of climate change on these compared to fungal communities. Final outcomes to cork oak forest sustainability will be dictated by the complex microbial network occurring between fungi and bacteria. More studies on specific microbial interactions will provide valuable information about which can be further explored for the protection of cork oak forests.

Availability of data and materials

Raw reads for each sample from both microbial communities were deposited in NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under BioProject number PRJNA744048.

References

  1. Jansson JK, Hofmockel KS. Soil microbiomes and climate change. Nat Rev Microbiol. 2020;18(1):35–46.

    Article  CAS  PubMed  Google Scholar 

  2. Meisner A, Jacquiod S, Snoek BL. Ten Hooven FC, van der Putten WH. Drought legacy effects on the composition of soil fungal and prokaryote communities. Front Microbiol. 2018;9(MAR):294.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Bardgett RD, Caruso T. Soil microbial community responses to climate extremes: resistance, resilience and transitions to alternative states. Philos Trans R Soc B Biol Sci. 2020;375:20190112.

    Article  CAS  Google Scholar 

  4. Elad Y, Pertot I. Climate change impacts on plant pathogens and plant diseases. J Crop Improv. 2014;28(1):99–139.

    Article  CAS  Google Scholar 

  5. Rolfe SA, Griffiths J, Ton J. Crying out for help with root exudates: adaptive mechanisms by which stressed plants assemble health-promoting soil microbiomes. Curr Opin Microbiol. 2019;49:73–82.

    Article  CAS  PubMed  Google Scholar 

  6. Vives-Peris V, de Ollas C, Gómez-Cadenas A, Pérez-Clemente RM. Root exudates: from plant to rhizosphere and beyond. Plant Cell Rep. 2020;39:3–17.

    Article  CAS  PubMed  Google Scholar 

  7. Fierer N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nat Rev Microbiol. 2017;15(10):579–90.

    Article  CAS  PubMed  Google Scholar 

  8. He D, Shen W, Eberwein J, Zhao Q, Ren L, Wu QL. Diversity and co-occurrence network of soil fungi are more responsive than those of bacteria to shifts in precipitation seasonality in a subtropical forest. Soil Biol Biochem. 2017;115:499–510.

    Article  CAS  Google Scholar 

  9. Preusser S, Poll C, Marhan S, Angst G, Mueller CW, Bachmann J, et al. Fungi and bacteria respond differently to changing environmental conditions within a soil profile. Soil Biol Biochem. 2019;137:107543.

    Article  CAS  Google Scholar 

  10. Sun Y, Chen HYH, Jin L, Wang C, Zhang R, Ruan H, et al. Drought stress induced increase of fungi:bacteria ratio in a poplar plantation. Catena. 2020;193:104607.

    Article  CAS  Google Scholar 

  11. Liu S, Wang H, Tian P, Yao X, Sun H, Wang Q, et al. Decoupled diversity patterns in bacteria and fungi across continental forest ecosystems. Soil Biol Biochem. 2020;144:107763.

    Article  CAS  Google Scholar 

  12. Bastías DA, Johnson LJ, Card SD. Symbiotic bacteria of plant-associated fungi: friends or foes? Curr Opin Plant Biol. 2020;56:1–8.

    Article  PubMed  CAS  Google Scholar 

  13. Martin FM, Uroz S, Barker DG. Ancestral alliances: plant mutualistic symbioses with fungi and bacteria. Science. 2017;356(819):eaad4501.

    Article  PubMed  CAS  Google Scholar 

  14. Deveau A, Bonito G, Uehling J, Paoletti M, Becker M, Bindschedler S, et al. Bacterial-fungal interactions: ecology, mechanisms and challenges. FEMS Microbiol Rev. 2018;42(3):335–52.

    Article  CAS  PubMed  Google Scholar 

  15. Lies A, Delteil A, Prin Y, Duponnois R. Using mycorrhiza helper microorganisms (MHM) to improve the mycorrhizal efficiency on plant growth. In: Role of Rhizospheric Microbes in Soil: Stress Management and Agricultural Sustainability. Singapore: Springer; 2018. p. 277–98.

  16. Reis F, Magalhães AP, Tavares RM, Baptista P, Lino-Neto T. Bacteria could help ectomycorrhizae establishment under climate variations. Mycorrhiza. 2021;1:3.

    Google Scholar 

  17. Jiang F, Zhang L, Zhou J, George TS, Feng G. Arbuscular mycorrhizal fungi enhance mineralisation of organic phosphorus by carrying bacteria along their extraradical hyphae. New Phytol. 2021;230(1):304–15.

    Article  CAS  PubMed  Google Scholar 

  18. Liu Y, Sun Q, Li J, Lian B. Bacterial diversity among the fruit bodies of ectomycorrhizal and saprophytic fungi and their corresponding hyphosphere soils. Sci Rep. 2018;8(1):1–10.

    Google Scholar 

  19. Trivedi P, Leach JE, Tringe SG, Sa T, Singh BK. Plant–microbiome interactions: from community assembly to plant health. Nat Rev Microbiol. 2020;18(11):607–21.

    Article  CAS  PubMed  Google Scholar 

  20. Lionello P, Scarascia L. The relation between climate change in the Mediterranean region and global warming. Reg Environ Chang. 2018;18(5):1481–93.

    Article  Google Scholar 

  21. Cramer W, Guiot J, Fader M, Garrabou J, Gattuso JP, Iglesias A, et al. Climate change and interconnected risks to sustainable development in the Mediterranean. Nat Clim Chang. 2018;8(11):972–80.

    Article  Google Scholar 

  22. Gauquelin T, Michon G, Joffre R, Duponnois R, Génin D, Fady B, et al. Mediterranean forests, land use and climate change: a social-ecological perspective. Reg Environ Chang. 2018;18(3):623–36.

    Article  Google Scholar 

  23. Moricca S, Linaldeddu BT, Ginetti B, Scanu B, Franceschini A, Ragazzi A. Endemic and emerging pathogens threatening cork oak trees: management options for conserving a unique forest ecosystem. Plant Dis. 2016;100(11):2184–93.

    Article  PubMed  Google Scholar 

  24. Bettenfeld P, Fontaine F, Trouvelot S, Fernandez O, Courty PE. Woody plant declines. What’s wrong with the microbiome? Trends Plant Sci. 2020;xx(xx):1–14.

    Google Scholar 

  25. Barrico L, Rodríguez-Echeverría S, Freitas H. Diversity of soil basidiomycete communities associated with Quercus suber L. in Portuguese montados. Eur J Soil Biol. 2010;46(5):280–7.

    Article  Google Scholar 

  26. Fumi MD, Mazzoleni V, Novelli E, Galli R, Busconi M, Blaghen M, et al. Pedologic characteristics and fungi community in unmanaged cork oak forest soil of two Mediterranean regions : Sardinia and Tunisia. Integr Prot Oak For IOBC WPRS Bull. 2014;101:31–8.

    Google Scholar 

  27. Maghnia FZ, Sanguin H, Abbas Y, Verdinelli M, Kerdouh B, El Ghachtouli N, et al. Impact of cork oak management on the ectomycorrhizal fungal diversity associated with Quercus suber in the Mâamora forest (Morocco). Comptes Rendus Biol. 2017;340(5):298–305.

    Article  Google Scholar 

  28. Reis F, Valdiviesso T, Varela C, Tavares RM, Baptista P, Lino-Neto T. Ectomycorrhizal fungal diversity and community structure associated with cork oak in different landscapes. Mycorrhiza. 2018;28(4):357–68.

    Article  PubMed  Google Scholar 

  29. Maghnia FZ, Abbas Y, Mahé F, Kerdouh B, Tournier E, Ouadji M, et al. Habitat- and soil-related drivers of the root-associated fungal community of Quercus suber in the northern Moroccan forest. PLoS One. 2017;12(11):1–17.

    Article  CAS  Google Scholar 

  30. Orgiazzi A, Lumini E, Nilsson RH, Girlanda M, Vizzini A, Bonfante P, et al. Unravelling soil fungal communities from different mediterranean land-use backgrounds. PLoS One. 2012;7(4):e34847.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Tardy V, Spor A, Mathieu O, Lévèque J, Terrat S, Plassart P, et al. Shifts in microbial diversity through land use intensity as drivers of carbon mineralization in soil. Soil Biol Biochem. 2015;90:204–13.

    Article  CAS  Google Scholar 

  32. Dubey A, Malla MA, Khan F, Chowdhary K, Yadav S, Kumar A, et al. Soil microbiome: a key player for conservation of soil health under changing climate. Biodivers Conserv. 2019;28(8–9):2405–29.

    Article  Google Scholar 

  33. Emberger L. Sur une formule climatique et ses applications en botanique. La Meteorol. 1932;92:423–32.

    Google Scholar 

  34. Emberger L. Une classification biogéographique des climats. In: Recueil des Travaux des Laboratoires de Botanique. Série Bota. Montpellier: Recueil des Travaux des Laboratoires de Botanique, Géologie et Zoologie de la Faculté des Sciences de L’Université de Montpellier; 1955. p. 3–43.

    Google Scholar 

  35. Haghighi AT, Zaki NA, Rossi PM, Noori R, Hekmatzadeh AA, Saremi H, et al. Unsustainability syndrome-from meteorological to agricultural drought in arid and semi-arid regions. Water (Switzerland). 2020;12(3):838.

    Google Scholar 

  36. Abatzoglou JT, Dobrowski SZ, Parks SA, Hegewisch KC. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015. Sci Data. 2018;5(1):1–12.

    Article  Google Scholar 

  37. QGIS Development Team. QGIS: a free and open source Geographic Information (version 13.14.16-Pi). 2020. http://www.qgis.org.

  38. Mitchell JI, Zuccaro A. Sequences, the environment and fungi. Mycologist. 2006;20(2):62–74.

    Article  Google Scholar 

  39. Chelius MKK, Triplett EWW. The diversity of archaea and bacteria in association with the roots of Zea mays L. Microb Ecol. 2001;41(3):252–63.

    Article  CAS  PubMed  Google Scholar 

  40. Tedersoo L, Bahram M, Põlme S, Kõljalg U, Yorou NS, Wijesundera R, et al. Global diversity and geography of soil fungi. Science. 2014;346(6213):1256688.

    Article  PubMed  CAS  Google Scholar 

  41. Bodenhausen N, Horton MW, Bergelson J. Bacterial communities associated with the leaves and the roots of Arabidopsis thaliana. PLoS One. 2013;8(2):e56329.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Illumina Inc. 16s metagenomic sequencing library preparation. 2013. https://support.illumina.com/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf.

  43. Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinformatics. 2011;27(6):863–4. Available from: http://prinseq.sourceforge.net/.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Joshi NA, Fass JN. Sickle: a sliding-window, adaptive, quality-based trimming tool for FastQ files (version 1.33). 2011. https://github.com/najoshi/sickle.

  45. Nurk S, Bankevich A, Antipov D, Gurevich A, Korobeynikov A, Lapidus A, et al. Assembling genomes and mini-metagenomes from highly chimeric reads. In: M D, Jiang R, Sun F, Zhang X, editors. Research in Computational Molecular Biology. vol 7821. Berlin: Springer; 2013. p. 158–70.

    Chapter  Google Scholar 

  46. Edgar RCR. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26(19):2460–1. Available from: http://www.drive5.com/usearch.

    Article  CAS  Google Scholar 

  47. Aronesty E. ea-utils: command-line tools for processing biological sequencing data. 2011. https://github.com/ExpressionAnalysis/ea-utils.

  48. Albanese D, Fontana P, De FC, Cavalieri D, Donati C. MICCA: a complete and accurate software for taxonomic profiling of metagenomic data. Sci Rep. 2015;5(1):1–7.

    Article  CAS  Google Scholar 

  49. Kõljalg U, Nilsson HR, Schigel D, Tedersoo L, Larsson K-H, May TW, et al. The taxon hypothesis paradigm—on the unambiguous detection and communication of taxa. Microorganisms. 2020;8(12):1910.

    Article  PubMed Central  CAS  Google Scholar 

  50. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41(Database issue):D590–6.

    CAS  PubMed  Google Scholar 

  51. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7(5):335–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. R core team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2019. Available from: https://www.r-project.org/.

    Google Scholar 

  53. Lahti L, Shetty S, Turaga N, Leung E, Gilmore R, Salojärvi J, et al. Tools for microbiome analysis in R. Microbiome package (version 1.17.42). 2017. http://microbiome.github.com/microbiome.

  54. Hughes JB, Bohannan BJM. Section 7 update: application of ecological diversity statistics in microbial ecology. In: Kowalchuk GA, de Bruijn FJ, Head IM, Akkermans AD, van Elsas JD, editors. Molecular Microbial Ecology Manual. 2nd ed. Dordrecht: Springer Netherlands; 2008. p. 3223–46.

    Chapter  Google Scholar 

  55. Magurran AE. Ecological diversity and its measurement. In: Ecological Diversity and Its Measurement. Dordrecht: Springer Netherlands; 1988.

    Chapter  Google Scholar 

  56. Hill TCJ, Walsh KA, Harris JA, Moffett BF. Using ecological diversity measures with bacterial communities. FEMS Microbiol Ecol. 2003;43(1):1–11.

    Article  CAS  PubMed  Google Scholar 

  57. Kassambara A. ggpubr: “ggplot2” Based Publication Ready Plots (version 0.4.0). 2020. https://cran.r-project.org/package=ggpubr.

  58. McMurdie PJ, Holmes S. Phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. Watson M, editor. PLoS One. 2013;8(4):e61217.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: community ecology package (version 2.5–7). 2020. https://CRAN.R-project.org/package=vegan.

  60. Kruskal JB, Wish M. In: Uslaner EM, editor. Multidimensional scaling. 07–001 ed. Beverly Hills and London: SAGE Publications; 1978.

    Chapter  Google Scholar 

  61. Henderson PA, Seaby RMH. Community analysis package. Lymington: Pisces Conservation Ltd; 2014.

    Google Scholar 

  62. Chambers JM. SoDA: functions and examples for “software for data analysis” (version 1.0-6.1). 2020. https://cran.r-project.org/package=SoDA.

  63. Macia-Vicente JG. comecol: Community ecology analysis (version 0.0.0.9000). 2020. https://github.com/jgmv/comecol.

  64. Warnes GR, Bolker B, Bonebakker L, Gentleman R, Huber W, Liaw A, et al. gplots: Various R programming tools for plotting data (version 3.1.1). 2020. https://cran.r-project.org/package=gplots.

  65. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12(6):R60.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Harrell Jr FE. Hmisc: Harrell Miscellaneous. 2021. https://cran.r-project.org/package=Hmisc.

  67. Bates D, Maechler M. Matrix: sparse and dense matrix classes and methods. 2019. https://cran.r-project.org/package=Matrix.

  68. Bastian M, Heymann S, Jacomy M. Gephi: an open source software for exploring and manipulating networks. In: Proceedings of the International AAAI Conference on Weblogs and Social Media. San Jose: AAAI Press; 2009;3(1):361–2.

  69. Nguyen NH, Song Z, Bates ST, Branco S, Tedersoo L, Menke J, et al. FUNGuild: an open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecol. 2016;20:241–8.

    Article  Google Scholar 

  70. Louca S, Parfrey LW, Doebeli M. Decoupling function and taxonomy in the global ocean microbiome. Science. 2016;353(6305):1272–7.

    Article  CAS  PubMed  Google Scholar 

  71. Wei T, Simko V. R package 'corrplot': visualization of a Correlation Matrix (version 0.84). 2017. https://github.com/taiyun/corrplot.

  72. Hofer U. The majority is uncultured. Nat Rev Microbiol. 2018;16:716–7.

    Article  CAS  PubMed  Google Scholar 

  73. Hurdeal VG, Gentekaki E, Hyde KD, Jeewon R. Where are the basal fungi? Current status on diversity, ecology, evolution, and taxonomy. Biologia (Bratisl). 2021;76(2):421–40.

    Article  Google Scholar 

  74. Carson JK, Gonzalez-Quiñones V, Murphy DV, Hinz C, Shaw JA, Gleeson DB. Low pore connectivity increases bacterial diversity in soil. Appl Environ Microbiol. 2010;76(12):3936–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Reis F, Soares-Castro P, Costa D, Tavares RM, Baptista P, Santos PM, et al. Climatic impacts on the bacterial community profiles of cork oak soils. Appl Soil Ecol. 2019;143(November 2018):89–97.

    Article  Google Scholar 

  76. Větrovský T, Kohout P, Kopecký M, Machac A, Man M, Bahnmann BD, et al. A meta-analysis of global fungal distribution reveals climate-driven patterns. Nat Commun. 2019;10(1):1–9.

    Article  CAS  Google Scholar 

  77. Talley SM, Coley PD, Kursar TA. The effects of weather on fungal abundance and richness among 25 communities in the intermountain west. BMC Ecol. 2002;2(1):1–11.

    Article  Google Scholar 

  78. Zhao Q, Jian S, Nunan N, Maestre FT, Tedersoo L, He J, et al. Altered precipitation seasonality impacts the dominant fungal but rare bacterial taxa in subtropical forest soils. Biol Fertil Soils. 2017;53(2):231–45.

    Article  Google Scholar 

  79. Reis F, Tavares RM, Baptista P, Lino-Neto T. Mycorrhization of Fagaceae forests within Mediterranean ecosystems. In: Mycorrhiza - function, diversity, state of the art. 4th ed. New York: Springer International Publishing; 2017. p. 75–97.

    Chapter  Google Scholar 

  80. Baptista P, Reis F, Pereira E, Tavares RM, Santos PM, Richard F, et al. Soil DNA pyrosequencing and fruitbody surveys reveal contrasting diversity for various fungal ecological guilds in chestnut orchards. Environ Microbiol Rep. 2015;7(6):946–54.

    Article  CAS  PubMed  Google Scholar 

  81. Wu BW, Gao C, Chen L, Buscot F, Goldmann K, Purahong W, et al. Host phylogeny is a major determinant of fagaceae-associated ectomycorrhizal fungal community assembly at a regional scale. Front Microbiol. 2018;9:2409.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Maghnia FZ, Abbas Y, Mahé F, Prin Y, El Ghachtouli N, Duponnois R, et al. The rhizosphere microbiome: a key component of sustainable cork oak forests in trouble. For Ecol Manag. 2019;434(December 2018):29–39.

    Article  Google Scholar 

  83. Habtewold JZ, Helgason BL, Yanni SF, Janzen HH, Ellert BH, Gregorich EG. Warming effects on the structure of bacterial and fungal communities in diverse soils. Appl Soil Ecol. 2021;163(August 2020):103973.

    Article  Google Scholar 

  84. Truong C, Gabbarini LA, Corrales A, Mujic AB, Escobar JM, Moretto A, et al. Ectomycorrhizal fungi and soil enzymes exhibit contrasting patterns along elevation gradients in southern Patagonia. New Phytol. 2019;222(4):1936–50.

    Article  CAS  PubMed  Google Scholar 

  85. Ren C, Liu W, Zhao F, Zhong Z, Deng J, Han X, et al. Soil bacterial and fungal diversity and compositions respond differently to forest development. Catena. 2019;181:104071.

    Article  CAS  Google Scholar 

  86. Suarez C, Ratering S, Schäfer J, Schnell S. Ancylobacter pratisalsi sp. nov. with plant growth promotion abilities from the rhizosphere of Plantago winteri Wirtg. Int J Syst Evol Microbiol. 2017;67(11):4500–6.

    Article  CAS  PubMed  Google Scholar 

  87. Lindström K, Mousavi SA. Effectiveness of nitrogen fixation in rhizobia. Microb Biotechnol. 2020;13(5):1314–35.

    Article  PubMed  CAS  Google Scholar 

  88. Dukunde A, Schneider D, Schmidt M, Veldkamp E, Daniel R. Tree species shape soil bacterial community structure and function in temperate deciduous forests. Front Microbiol. 2019;10:1519.

    Article  PubMed  PubMed Central  Google Scholar 

  89. de Vries FT, Griffiths RI, Bailey M, Craig H, Girlanda M, Gweon HS, et al. Soil bacterial networks are less stable under drought than fungal networks. Nat Commun. 2018;9(1):3033.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  90. Shen C, Gunina A, Luo Y, Wang J, He JZ, Kuzyakov Y, et al. Contrasting patterns and drivers of soil bacterial and fungal diversity across a mountain gradient. Environ Microbiol. 2020;22(8):3287–301.

    Article  PubMed  Google Scholar 

  91. Yu WC, Zhou X, Guo D, Hua ZJ, Yan L, Zhong FG, et al. Soil pH is the primary factor driving the distribution and function of microorganisms in farmland soils in northeastern China. Ann Microbiol. 2019;69(13):1461–73.

    Article  CAS  Google Scholar 

  92. Xu X, Chen C, Zhang Z, Sun Z, Chen Y, Jiang J, et al. The influence of environmental factors on communities of arbuscular mycorrhizal fungi associated with Chenopodium ambrosioides revealed by MiSeq sequencing investigation. Sci Rep. 2017;7(1):1–11.

    CAS  Google Scholar 

  93. Seaton FM, George PBL, Lebron I, Jones DL, Creer S, Robinson DA. Soil textural heterogeneity impacts bacterial but not fungal diversity. Soil Biol Biochem. 2020;144:107766.

    Article  CAS  Google Scholar 

  94. Vera A, Moreno JL, Siles JA, López-Mondejar R, Zhou Y, Li Y, et al. Interactive impacts of boron and organic amendments in plant-soil microbial relationships. J Hazard Mater. 2021;408:124939.

    Article  CAS  PubMed  Google Scholar 

  95. Vera A, Moreno JL, García C, Morais D, Bastida F. Boron in soil: the impacts on the biomass, composition and activity of the soil microbial community. Sci Total Environ. 2019;685:564–73.

    Article  CAS  PubMed  Google Scholar 

  96. Wu L, Ma H, Zhao Q, Zhang S, Wei W, Ding X. Changes in soil bacterial community and enzyme activity under five years straw returning in paddy soil. Eur J Soil Biol. 2020;100:103215.

    Article  CAS  Google Scholar 

  97. Xu Z, Zhang T, Wang S, Wang Z. Soil pH and C/N ratio determines spatial variations in soil microbial communities and enzymatic activities of the agricultural ecosystems in Northeast China: Jilin Province case. Appl Soil Ecol. 2020;155:103629.

    Article  Google Scholar 

  98. Luo Y, Iqbal A, He L, Zhao Q, Wei S, Ali I, et al. Long-term no-tillage and straw retention management enhances soil bacterial community diversity and soil properties in southern China. Agronomy. 2020;10(9):1233.

    Article  CAS  Google Scholar 

  99. Soleimani A, Hosseini SM, Massah Bavani AR, Jafari M, Francaviglia R. Influence of land use and land cover change on soil organic carbon and microbial activity in the forests of northern Iran. Catena. 2019;177:227–37.

    Article  CAS  Google Scholar 

  100. Odriozola I, Abrego N, Tláskal V, Zrůstová P, Morais D, Větrovský T, et al. Fungal communities are important determinants of bacterial community composition in deadwood. Shade A, editor. mSystems. 2021;6(1):e01017–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Collado E, Bonet JA, Camarero JJ, Egli S, Peter M, Salo K, et al. Mushroom productivity trends in relation to tree growth and climate across different European forest biomes. Sci Total Environ. 2019;689:602–15.

    Article  CAS  PubMed  Google Scholar 

  102. Müller K, Kubsch N, Marhan S, Mayer-Gruner P, Nassal P, Schneider D, et al. Saprotrophic and ectomycorrhizal fungi contribute differentially to organic P mobilization in beech-dominated forest ecosystems. Front For Glob Chang. 2020;3:47.

    Article  Google Scholar 

  103. Maillard F, Kennedy PG, Adamczyk B, Heinonsalo J, Buée M. Root presence modifies the long-term decomposition dynamics of fungal necromass and the associated microbial communities in a boreal forest. Mol Ecol. 2021;30(8):1921–35.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

The authors thank Maria Carolina Varela for the help in identification of disease severity levels of trees and Pinto Moreira (HC forests), Eng. Carlos Barroso Dias (PNPG), Eng. Luís Dias (Herdade das Barradas da Serra—GR) and José Lino Neto (GV forest) for providing access to cork oak forests.

Funding

This work was supported by FEDER funds through COMPETE (Programa Operacional Factores de Competitividade) and by national funds by FCT (Fundação para a Ciência e a Tecnologia) in the framework of the project POCI-01-0145-FEDER-028635, as well as for financial support to BioISI (UIDB/04046/2020), CBMA (UIDB/04050/2020) and CIMO (UID/AGR/00690/2020) by national funds by FCT/MCTES/PIDDAC. D. Costa thanks FCT for PhD grant SFRH/BD/120516/2016 and COVID/BD/151779/2021.

Author information

Authors and Affiliations

Authors

Contributions

D.C. contributed to study design and collection of samples, performed the experiments, analyzed data, and contributed to manuscript writing. R.M.T. contributed to study design and manuscript editing. P.B. contributed to study design and manuscript editing. T.L.N. contributed to study design, collection of samples and manuscript writing. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Teresa Lino-Neto.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

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://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Costa, D., Tavares, R.M., Baptista, P. et al. The influence of bioclimate on soil microbial communities of cork oak. BMC Microbiol 22, 163 (2022). https://doi.org/10.1186/s12866-022-02574-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12866-022-02574-2

Keywords