Study area
The Hengduan Mountains region (HDMR) – located at the eastern edge of the Qinghai-Tibetan plateau, China – is a core region of the Sino-Himalayan Floristic Region and one of Earth's 34 biodiversity hotspots, where six parallel mountain chains and rivers stretch from north to south and the traffic from west to east is blocked (Fig. S1). HDMR has a total area of 500,000 km2 and includes the eastern Tibet Autonomous Region, the western Sichuan Province, and the northwestern Yunnan Province in which there are around 99 counties. The generalized geographic range of HDMR is between 22 and 33° N latitude and between 97 and 103° E longitude and starts from the Bayankala Mountains in the north to the Gaoligong Mountains in the south and from the Qionglai Mountains in the east to the Boshula Mountain ridge in the west, the latter representing the steering extension of the Nyainqentanglha and Tanggula mountains (Gan, 2007). The elevation ranges from 386 to 7143 m and is higher in the north-western area than the south-eastern area. The topography declines from north-west to south-east. As HDRM is located on a plateau covered by both a temperate and a subtropical climate zone, the climate differs distinctly due to the complex mountainous topography (Li et al., 2011; Zhu et al., 2012). Thus, with abundant precipitation in the accumulation area and low temperature in the higher altitude area, HDMR is the southern- and eastern-most glacial region in China (Li et al., 2011; Li et al., 2010). HDMR is also located in a typical monsoon climate zone that is controlled by the South Asia and the East Asia monsoons (Li & Su, 1996). In addition, HDMR covers both sides of the ‘Hu line’ which was proposed by the Chinese geographer Hu Huanyong and today is regarded as the Chinese population density distribution line (Fig. S1). At the north-western side of the Hu line, population density is sparse and economic development relatively slow, while the opposite is true for the south-eastern side.
Field sampling
Field surveys of macrophyte communities were conducted during the similar season (summer: June to August) 2014–2018. A total of 802 local macrophyte communities were sampled across 48 of 99 county in HDMR and the abundance of macrophytes was estimated using 3–5 randomly positioned quadrats (0.25 m2) at each site. Details on macrophyte sampling methods are given in Fu et al. (2017). All species occurring in the plots were identified and recorded. At each site, we measured sampling depth and water pH and recorded the geographic coordinates (i.e., latitude and longitude) and altitude using a portable GPS.
The HDMR consists of 99 county-level administrative division, and these counties are markedly disconnected from each other due to typical geographic barriers (big mountains, > 4000 m a.s.l). In this study, a suite of local communities (N ≥ 15, sites within or near a river or a lake, temporarily connected ponds) originating from a single county were defined a metacommunity. In total, the 802 local communities consisted of 48 metacommunities, each of which contained 15–77 local communities (Fig. S1).
Elements of metacommunity structure
To test which metacommunity structure best fitted the macrophyte species distribution, we used the EMS methodology described by Leibold & Mikkelson (2002) and Presley et al. (2010). Following the ‘range perspective’, metacommunity structure was assessed by hierarchically evaluating the coherence, turnover, and boundary clumping elements of a site using a species incidence matrix, which was ordinated using reciprocal averaging (i.e., correspondence analysis).
First, coherence was calculated as the number of embedded absences in an ordinated incidence matrix. Then, the significance of coherence was tested by comparing the observed values to a null distribution of 999 simulated matrices. A significant negative coherence (i.e., more absences than expected by chance) indicates the existence of a checkerboard metacommunity type with strong interspecific competition, while non-significant coherence is related to randomness of species distribution with regard to gradient. A significant positive coherence (i.e., less absence than expected by chance) suggests that a common gradient shapes the species distribution within a metacommunity (Leibold & Mikkelson, 2002), and the metacommunity type was further specified by assessing the turnover and boundary clumping patterns.
Second, turnover was calculated as the number of times that one species replaces (Rep) another between two sites in an ordinated incidence matrix. The significance of turnover was then tested by comparing the observed number of replacements to the ones calculated from the null distribution. A significantly lower turnover than expected indicates a nestedness pattern. A significantly higher turnover than expected suggests evenly spaced Gleasonian or Clementsian patterns that were identified with further analysis of boundary clumping.
Third, boundary clumping was evaluated as the Morisita’s Index (MI, (Morisita, 1971)), which was compared to a null expectation of 1 using an χ2 test. A randomly distributed range boundary (i.e., Gleasonian) is expected when MI does not differ from 1, and a clumped range boundary (i.e., Clementsian) is expected when MI is significantly > 1, and a hyperdispersed range boundary (i.e., evenly spaced metacommunity type) occurs when MI is significantly < 1.
In addition, quasi-structures of a metacommunity appear when significant positive coherence and nonsignificant turnover exist, and these are relatively weaker than the standardized ones. Thus, non-significant negative turnover suggests quasi-nestedness, and non-significant positive turnover suggests quasi-Gleasonian, quasi-Clementsian, or quasi-evenly spaced gradients, which can be separated by boundary clumping (Presley et al., 2010).
To detect the significance of coherence and turnover, we applied a fixed-proportional null model (r1) to create 999 random matrices with fixed species richness for each site (rows). We calculated the Z-score for coherence and the turnover for each metacommunity as follows:
Z-score = (Iobs – Inull) / Isd,
where Iobs is the observed index values for coherence (Abs) and turnover (Rep), Inull is the mean index values calculated from 999 null models, and Isd is the standard deviation of the index values calculated from 999 null models. We checked symmetry of null-distribution before calculating the SES value and used log-transformation Z-score of the test statistic when the null distribution is skewed (Botta-Dukát, 2016). Z-scores allow comparisons among metacommunities and can thus subsequently be used in comparative analyses. A Z-score higher than 1.97 and a Z-score lower than − 1.97 indicate significance at 0.05 level. EMS analysis was conducted using functions of the “metacom” package.
Statistical analysis
Local variables (e.g., sampling depth and water pH) and geographic variables (e.g., altitude and latitude) were calculated as polygon centroids across sites within each metacommunity. We calculated spatial extent within each metacommunity using the Euclidean distance between sites. In addition, we classified HDMR into two parts (north-west and south-east) according to the “Hu Line” as mentioned above, this line being considered a regional variable (hereafter called HL). All variables were log-/square root-transformed prior to the analysis if the distribution was not normal. We tested the correlations among the six variables using the Spearman method.
We were also interested in testing how alpha, beta, and gamma aspects of biodiversity measures were related to the metacommunity structure. Alpha diversity was calculated as the mean species richness across the sites within a metacommunity and gamma diversity as the sum of species occurring in a metacommunity. Total beta diversity (Sørensen coefficient) and its two components – turnover (Simpson coefficient) and nestedness (nestedness coefficient) – were calculated using the functions ‘beta.multi’ for multiple-site indices in the R package ‘betapart’ (Baselga & Orme, 2012).
In order to determine which variables contributed to explaining the geographic patterns of macrophyte metacommunities in HDMR, we used two complementary approaches. For the three elements of metacommunity structure, we first used generalized linear model (GLM) to assess the effects of environmental variables and diversity measures (predictors), applying separately the Z-scores of coherence, the Z-scores of turnover, or the index of boundary clumping as response variables. The variance inflation factors (VIF) of all tested variables were < 6.1, indicating that there was no problem of collinearity among the predictor variables. Next, we applied the model selection approaches using the second-order Akaike’s information criterion (delta AIC < 2) to select the best models with the most important explanatory variables for the beta diversity metrics. The sum of Akaike weights including all models was calculated to estimate the relative importance of explanatory variables. The model selection and model averaging were conducted using functions of the “MuMIn” package (Bartoń, 2018). We also applied the above statistical approaches to assess the relative importance of environmental variables (predictors) for determining the variations of the alpha, beta, and gamma aspects of biodiversity measures.
To distinguish metacommunity types, we used linear discriminant function analysis (DFA) to evaluate how well the six environmental variables and the three diversity measures maximized the differences in the structure of the 48 macrophyte metacommunities using all the detected categorical “metacommunity types” as response variable. DFA was conducted using the function “lda” in the R package MASS (Ripley et al., 2013). We also used stepwise selection of predictor variables to see which predictors were most important in separating the metacommunity types using the function “greedy.wilks” in the R package klaR (Weihs et al., 2005).
All statistical tests were performed using R version 3.51 software (Core, 2013). The figures were constructed using the R package “ggplot2” (Wickham, 2009).