Plant fungal diseases have a great influence on both photosynthesis and ecosystem function. However, how the elevation gradient, which is one of the most biogeographic factors, affects diseases is scarce. Here, we combined a field survey and a meta-analysis to test how elevation affect foliar fungal diseases and soil fungal pathogens through different paths. We arranged 30 plots along 3200 m ~ 4000 m in a Qinghai-Tibetan alpine meadow and collected the data of foliar fungal diseases, plant composition and soil properties to study how environment-mediated (through changes in temperature and humidity), plant community-mediated (through changes in plant biomass, richness, evenness, phylogenetic structure and community composition) and soil mediated (through changes in soil properties) effects of elevation on foliar fungal diseases and soil fungal pathogens. Based on linear models, we found that elevation decreased soil fungal pathogen richness rather than community pathogen load of foliar diseases. Specifically, a combination of community proneness and Pielou’s evenness index was the best model in predicting pathogen load. The structural equation model further confirmed that although elevation significantly changed both the plant community indices and soil properties, elevation mainly drove pathogen load via plant community-mediated effects, but decreased soil fungal pathogen richness through temperature. A systematic meta-analysis composed of 48 studies from 31 literatures confirmed our main conclusions that elevation did not significantly foliar fungal diseases, but decreased soil fungal pathogen richness significantly, indicating contrast effects of elevation in driving above- and belowground plant pathogens. Hence, we distinguished the different mechanisms for different parts of the plant pathogens in one system, and our study will improve the predictability of plant diseases, especially under the background of global climate change.

Field survey along elevation

Study site and plots establishment

Our survey was conducted along an elevation gradient ranging from 3200 m to 4000 m on the south slope of Qilian Mountains in Menyuan County, located in the northeastern Qinghai-Tibetan Plateau. The climate is continental monsoon climate, with a 6 months growing season from mid-April to mid-October (Ma et al., 2017). The vegetation is alpine meadow, generally dominated by some genera of perennial herbaceous (e.g. Gentiana, Kobresia, Poa and Saussurea), with the species composition shift along elevation. We established thirty 0.5 × 0.5 m plots (6 replications for each of five elevations) along the elevation gradient: 3200 m (37°36'39'' N, 101°18'16'' E), 3400 m (37°39'58'' N, 101°20'20'' E), 3600 m (37°41'47'' N, 101°21'34'' E), 3800 m (37°42'13'' N, 101°22'14'' E) and 4000 m (37°42'29'' N, 101°22'27'' E). The plots were randomly selected with at least 10 m buffer zone between two adjacent plots at each elevation.

Plant and soil sampling and climate monitoring

In early August 2020, we harvested all the plant aboveground parts at ground level and sorted them into plant species at plot-level. Then we dried them at 65℃ for 48 hours to constant mass and weighed to 0.01 g, and sumed the dry biomass from a plot as the aboveground biomass (hereafter as ‘AB’, see Table S1 for abbreviations). We also collected four soil cores (5 cm in diameter and 10 cm in depth) from each plot and pooled them as one sample. All fine roots were collected from each sample and were flushed by water, then we dried them for plant belowground biomass (BB). Soil moisture content (W; %) was measured gravimetrically after 5 h of desiccation at 120℃. A pH analyzer and a conductivity analyzer were used to measure the soil pH (pH) and soil conductivity (C; ms/s), respectively. 5 grams of fresh soil were extracted with 50 ml 0.2 M KCl for 1 h at 60 rev s^-1 using a shaker, then nitrate-nitrogen (NO₃^-; mg/kg) and ammonium-nitrogen (NH₄⁺; mg/kg) were measured using an auto-analyzer (AA3, Bran-Luebbe, Germany).

We recorded foliar fungal disease severity following the methods provided in Liu et al. (2017). In brief, we recorded disease severity (i.e. % leaf area coverd by fungal lesion; V_i) from five leaves randomlty selected from five individuals for each plant species in each plot. We recorded all available leaves for species with less than 25 leaves. Pathogen identification mainly followed identification manuals including Fungal Identification Manual (Wei, 1979), and also previous studies in this area (Zhang, 2009; Liu et al., 2019; Liu et al., 2020b).

We placed a Temperature-Humidity Recorder Cos-03-0 (Renke Control Technology Co., Ltd., Jinan, Shandong, China) at each elevation, which continuously monitored the ambient temperature and humidity every 60 seconds for nonstop 72 hours in early August. Then we calculated the mean daily temperature (MDT) and mean daily humidity (MDH) for each elevation.

DNA extraction and PCR amplification

Soil total DNA was extracted by using a TIANGEN Magnetic Soil And Stool DNA Kit (TIANGEN Biotech Co., Ltd., Beijing, China) following the manufacturer’s protocol. The concentration and purity of the DNA were accessd by agarose gel electrophoresis. The primers internal transcribed spacer 1 (ITS1) region was amplified with the forward primer ITS1F (5’-GGAAGTAAAAGTCGTAACAAGG-3’) and the reverse primer ITS1R (5’-GCTGCGTTCTTCATCGATGC-3’) using polymerase chain reactions (PCRs) (Gardes & Bruns, 1993). PCRs were performed in a 30 μL mixture composed of 15 µL of Phusion Master Mix (2 ×), 10 µL of DNA template, 3 µL of primer, and 2 µL of dd H₂O. The PCR reactions were conducted as follows: a denaturation at 98 ℃ for 1 min, followed by 30 cycles of 30 s at 98 ℃, 30 s at 55 ℃ for annealing and 30 s at 72 ℃ for elongation, finally 5 mins at 72 ℃ for extention and reaction termination.

High-throughput sequencing and processing of the sequencing data

The library was constructed by a TruSeq DNA PCR-Free Library Preparation Kit (Illumina, San Diego, CA), then we used the Qubit@ 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA) and Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA) to assesse the quality of the library. The Illumina NovaSeq 6000 platform (Illumina, San Diego, CA) was operated to quantified and sequenced the constructed library (Novogene Co., Ltd., Beijing, China). Then non-clustering direct denoising to generate OTU based on unoise3 (UNBIAS) (https://www.drive5.com/usearch/manual/unoise_algo.html), and rare sequences (reads < 8) were removed to avoid the possible spurious reads generated by sequencing errors. UNITE 2020 (http://unite.zbi.ee; Kõljalg et al., 2005) was used to cluster the sequences into operational taxonomic units (OTUs).

A fungal species can be defined as a plant fungal pathogen when it belonged to a genus with the majority of their species can induce plant disease symptoms (Liang et al., 2016), referenced to the published data (Tedersoo et al., 2014) and FUNGuide database (https:// github.com/UMNFuN/FUNGuild; Nguyen et al., 2015). Then we calculated the accumulated OTU number of soil fungal pathogens (sfpOTUs), and also their relative abundance (sfpRA; sequences of soil fungal pathogens divided by the total number of sequences) for each sample respectively. Soil fungal pathogens relative abundance (sfpRA) was log-transformed to achieve the normality of residuals in the following analysis.

Community pathogen load

Community pathogen load (PL) was defined as (Mitchell et al., 2002):

where S is the total number of plant species in a certain plot, b_i is the aboveground biomass of plant species i.

We defined a ‘disease proneness index’ (hereinafter ‘P_i’) for each species as the average severity index (V_i) across 30 plots of plant species i. We then calculated a ‘community proneness index’ (hereinafter ‘p’) for each plot by calculating a plant aboveground biomass-weighted average of the P_i for each plot (Liu et al., 2017):

where p is the expected community pathogen load base on constituent host plants (Liu et al., 2017). Community pathogen load (PL) and disease proneness index (p) were log-transformed to achieve the normality of residuals in the following analysis.

Analysis

(1) Soil properties We used rda function in vegan package (Oksanen et al., 2020) to conduct principal component analysis (PCA) to summarize the soil properties (W, pH, C, NO₃^-, NH₄⁺), with the first principal component (Soil PCA1) explained 41.63% of the total variance (Table S2). Then we calculated the Spearman rank-order correlation between Soil PCA1 and each of soil properties (W, pH, C, NO₃^-, NH₄⁺) by using cor.test function. Specifically, the Soil PCA1 is negatively related to C, pH and NH₄⁺, while positively related to W and NO₃^-.

(2) Plant community-level indices For each plot, we calculated the plant species richness (SR) and Pielou's evenness index (Ev) using diversity function in vegan package. Then we built a plant phylogenetic tree based on rbcL and matK genes, following the methods provided in Liu et al. (2015). In brief, MEGA X (Kumar et al., 2018) was used to conduct multiple sequence alignment by employing ClustalW method (Larkin et al., 2007). Then we concatenated the aligned two genes sequences into a matrix. For plant species that missing rbcL or matK sequences, we used their congeneric representatives in NCBI (National Center for Biotechnology Information). We then selected the top-ranked maximum-likelihood model of nucleotide substitution using Akaike’s information criterion (AIC) by jModelTest v2.1.7 (Darriba et al., 2012), and built maximum-likelihood phylogenies using PhyML v3.1 (Guindon et al., 2010) with the aligned and concatenated sequences and also the optimal nucleotide substitution model (GTR + G). A 500 bootstrap replicates were run to evaluate the nodal support on maximum-likelihood phylogenies. We calculated several phylogenetic indices by picante package (Kembel et al., 2010), including standardized effect size of Faith's phylogenetic diversity, mean pairwise distance (weighted by aboveground biomass), mean nearest taxon distance (weighted by aboveground biomass; MNTD). Among these indices, MNTD was the best predictor for community pathogen load (PL), hence, we only used MNTD in the following analysis, given the high collinearity among these phylogenetic indices.

(3) Factors affecting foliar fungal diseases and soil pathogens At plant population-level, we evaluated how elevation (hereafter Ele) affects V_i of each plant species based on a series of linear models, and calculated the corresponding F- and P-values.

We set the Ele as independent variable, PL and soil pathogen indices (sfpOTUs, sfpRA) as response variables in a series of linear models to test the direct effect of Ele. Due to the collinearity among various community-level indices (SR, Ev, p, MNTD, AB, BB), environmental factors (MDT, MDH) and soil properties (Soil PCA1) (Fig. S1), we set them as independent variables in a series of linear models to test their effects on PL and various soil pathogen indices. Then we compared these models with the null model (i.e. the intercept-only model) based on Akaike's information criterion corrected for small sample sizes (AIC_c). We also calculated the log-likelihood (LL) and AIC_c based parameters: change in AIC_c relative to the top-ranked model (DAIC_c), AIC_c weight (wAIC_c) and the percent deviance explained (De) (Burnham et al., 2011), to estimate their possibilities to be the best predictor in predicting PL and various soil pathogen indices.

We plotted the correlation matrix for all biotic or abiotic variables (Ele, SR, Ev, p, MNTD, AB, BB, MDT, MDH, Soil PCA1), then we calculated the Pearson’s correlation between these variables and PL/various soil pathogen indices (Fig. S1a). We also conducted Mantel test based “bray-curtis” distance between sfpOTUs and biotic or abiotic variables using the ggcor package (Huang et al., 2020) (Fig. S1b). The Permutational multivariate analysis of variance (PERMANOVA) was conducted to test the compositional difference of soil fungal pathogens along the elevation gradient. We then established a piecewise structural equation model (piecewise SEM; Lefcheck, 2016) to test how elevation affects PL and soil pathogen indices through community-mediated, environment-mediated and soil-mediated effects (Fig. S7). Our piecewise SEM comprised a series of linear models. To simplify the paths, we only retained Ev to character plant diversity, because Ev was the best predictor for PL. And we retained sfpOTUs for soil pathogens, given it was sensitive to environmental factors. We calculated the standardized path coefficients (scaled by their mean and standard deviation) and corresponding significance (P values) for each path of the final models. The Fisher’s C test and AIC were used to test the goodness-of-fit of piecewise SEM by using piecewiseSEM package (Lefcheck, 2016). We further identified the effects of each index in piecewise SEM (Ev, p, Soil PCA1, MDT) on PL and sfpOTUs after accounting for the potential impact from others by residual plots.

Meta-analysis

Data collection

To test the universality of our main results of the field survey, we conducted a systematic literature search in ISI Web of Science and China National Knowledge Infrastructure (www.cnki.net). We did the search for foliar fungal diseases [(plant disease* OR pathogen* OR infect* OR epidemic*) AND (inciden* OR prevalen* OR load* OR severity OR occur* OR abundance) AND (elevation* OR altitud*) AND (tree* OR forest* OR grass* OR shrub*)] and soil plant pathogens [(fungal OTU* OR fungi OTU* OR fung* operational taxonomic unit OR fung* abundance) AND (elevation* OR altitud*) AND (soil OR belowground OR underground)] respectively. We initially obtained a set of 2680 publications, with 2168 for foliar fungal diseases and 512 for soil plant pathogens. We included only those studies that: (i) focusing on the relationship between elevation and foliar fungal diseases/soil plant pathogens in nonagricultural ecosystems; (ii) reporting more than three sample sizes. Finally, we indentified 48 studies from 31 literatures met our criteria (Table S8).

We collected the OTU table for studies focued on soil plant pathogens, identified the plant pathogens according to abovementioned methods, and calculated sfpOTUs for each of study. We extracted the sample sizes and Pearson’s correlation coefficients (r) from the main taxt, tables, figures (by WebPlotDigitizer Version 4.4; Rohatgi, 2020) or the raw data. We also recorded the background information of lacation, mean annual temperature, mean annual precipitation and the elevation span of sampling.

Effect sizes

We calculated the Fisher’s z-transformation of Pearson’s correlation coefficients (r) as effect size (Rosenberg et al., 2013):

the corresponding variance for each Z was estimated as:

where n is the sample size. Positive values of Z indicate elevation increases plant diseases/pathogens, while negative for decrease.

Analysis

We used the rma.mv function in the metafor package (Viechtbauer, 2010) to calculated the mean effect size of elevation on foliar fungal diseases/sfpOTUs, with ‘study’ nested in ‘paper’ as random structure (Nakagawa et al., 2017). The effect size (Z) was considered to be significant when the 95% confidence interval of the mean did not include zero (Lajeunesse, 2013). We tested the overall effect of elevation on foliar fungal diseases and sfpOTUs respectively. Then introduced mean annual temperature, mean annual precipitation, latitude and elevation span of sampling to test the context dependence of effect size (Z). The amount of heterogeneity explained by each variable was estimated by Q_m statistic and its corresponding P value (Viechtbauer, 2010). For assessing the potential publication bias, we conducted Kendall’s rank test for funnel plot asymmetry (Borenstein et al., 2009), and also did the meta-regression between effect size (Z) and literatures’ published year/journal impact factor. All statistical analyses were conducted using R v4.1.1 (R Development Core Team, 2021).

(1) HaibeiData.xlsx

This .xlsx file contains the raw data analyzed in the "Field survey" (which includes the data of plant community, Soil PCA1, environmental foctors and disease indices for each plot).

(2)MetaData.xlsx

This .xlsx file contains the raw data analyzed in the "Meta-analysis"(which includes the citation information, geographic information (location), research object information and statistical indices for each study).