**1. Introduction**

The forest–steppe ecotone is located at the driest edges of forest distribution and is continually threatened by the increasing temperature and drought [1–4]. Forest patches and grassland patches

distribute in a mosaic pattern in the forest–steppe ecotone [4]. To answer whether the forests are stable in the forest–steppe ecotone, it is helpful to explore what interspecific associations exist among tree species and what physiological mechanisms are driving them, as the interspecific associations of co-occurring trees will determine the forest composition variations [5].

Interspecific association—referring to the relationships and interactions among co-occurring species [6]—has been reported to be driven by many factors, including the forest succession process [5], the heterogeneity of microenvironment [7], etc. In the forest–steppe ecotone—the driest edges of forest distribution—water deficit is the most dominant environment stress affecting forest dynamics and functions [2,8], instead of reported interspecific association drivers. Moreover, water deficit here has been predicted to be exacerbated in the future [2,8]. For this reason, it is crucial to explore the importance and role of water deficit in shaping the vegetation community in the forest–steppe ecotone.

Facing the aggravated water deficit, plant hydraulic strategies balancing hydraulic safety and water-transport efficiency are believed to determine tree drought responses [9–11]. The water-transport efficiency of trees has direct impacts on the carbon uptake and productivity of trees, while xylem resistance to cavitation determines the safety of individuals when extreme drought events occur [11–13]. The safety and efficiency of trees generally occur as a tradeoff. A larger conduit diameter provides higher hydraulic efficiency as well as higher vulnerability to drought-induced embolism in the xylem, and vice versa [14,15]. The hydraulic safety–efficiency combination of trees can provide a way to map hydraulic strategy of each tree species, to describe the differences among the hydraulic strategies of co-occurring species and to judge the vulnerability of trees in the forest–steppe ecotone under ongoing climate change [16,17].

If water deficit plays a considerable role in forest interspecific associations, divergent hydraulic strategies should differentiate the water-use strategies of species and reduce interspecific water competitions [18], thus benefitting species co-occurrence. Further, we can expect the differences of species hydraulic strategies to be an effective predictor of interspecific relationships, i.e., the hydraulic strategies of co-existent species to be divergent, while the hydraulic strategies of mutually exclusive species to be convergent. Meanwhile, the species with both higher hydraulic safety and efficiency should have higher potential to be dominant in the community. Forest–steppe ecotone has all the ten genera of tree species distributing in semi-humid forests in northern China, even with a mean annual precipitation of about 300 mm [19]. If the diversity of tree hydraulic strategies can provide partitioning of the water resource, it may be the reason for forest–steppe ecotone to maintain the high tree species diversity. In general, we developed two hypotheses in this article, 1) there is a positive correlation between interspecific association and distance in species hydraulic strategy, and 2) species with both high hydraulic efficiency and safety may play a dominant role in the forest, or have experienced a rapid population development over the years.

To verify the hypotheses, we carried out a two-year vegetation survey and a series of hydraulic measurements in a permanent plot in the Saihan Wula National Nature Reserve, Inner Mongolia, China. In this work, hydraulic efficiency is represented by leaf-specific hydraulic conductivity (Kl), while hydraulic safety is represented by a series of important thresholds on the stem-vulnerability curves, stem water potential at a 50% loss of stem conductivity (P50) especially. The Kl and P50 are combined to show the hydraulic strategy of each species.

#### **2. Materials and Methods**

#### *2.1. Site Description and Sampling Design*

Field work was conducted within a 6-ha permanent plot in a temperate deciduous broadleaved forest in the Saihan Wula National Nature Reserve in the southern Da Hinggan Mountains, Inner Mongolia, China (Figure 1a,b, 44◦12 N, 118◦44 E). The plot is located in semiarid forest–steppe ecotone, with a mean annual temperature of 2 ◦C and a mean annual precipitation of 358 mm. The wet season is from May to September. The elevation of the plot is approximately 1175–1355 m ASL [20]. The studied forest is located on the north side of the mountain, while the steppe lies on the south (Figure 1d). Severe drought-induced forest mortality occurred during 2011–2012, followed by a huge amount of regeneration. Tree density increased from 1825.17 ha−<sup>1</sup> in 2012 to 3812.17 ha−<sup>1</sup> in 2015. No grazing, agricultural abandonment, forest managemen<sup>t</sup> or fire interference have been found inside or around the plot since 1997.

**Figure 1.** Study area location and site characteristics. (**a**) Location of the study area in China; (**b**) landscape of semiarid forest–steppe ecotone in which the study area lies; (**c**) sketch of the plot (shown in green polygon) and quadrats containing sampled trees for hydraulic architecture measurements (shown as yellow squares). Each quadrat has an area of 100 m2; (**d**) landscape photography of the study area showing the mosaic distribution of forest (dark green patches) and grassland (light green patches); (**e**) amount of each tree species in 2012 (shown in blue) and 2015 (shown in orange); (**f**) schematic picture showing the forest structure within the plot.

Two comprehensive plot surveys were carried out in the whole permanent plot in 2012 and 2015, respectively. The permanent plot was divided into 608 contiguous 10 m × 10 m quadrats. Basic information of all trees with diameter at breast height >1 cm was carefully recorded, including species, tree height, diameter at breast height and the located quadrat identification number. Only main species were taken into further analysis, with the proportion in the total tree number over 1%, which are *Populus davidiana* Dode, *Betula platyphylla* Suk., *Betula dahurica* Pall., *Quercus mongolica* Fisch. ex Ledeb and *Tilia mongolica* Maxim. (Figure 1e,f). *P. davidiana*, *B. platyphylla* and *B. dahurica* are early successional species, while *Q. mongolica* and *T. mongolica* are late successional species. *P. davidiana* was the most dominant tree species in the forests, accounting for the largest amount of the tree mortality and regeneration during and after the extensive drought events in 2011–2012 (Figure 1e), with a mortality ratio of over 60% [21].

Hydraulic traits were measured for five trees per species at different elevations (Figure 1c). The maximal vessel lengths of each species are: *P. davidiana*: 11 cm, *B. platyphylla*: 14 cm, *B. dahurica*: 16 cm, *Q. mongolica*: 86 cm and *T. mongolica*: 20 cm.

#### *2.2. Interspecific Association Quantification*

The interspecific association of five species were quantified with the plot survey data. The variance ratio (VR) was used to test the overall association among species. If VR = 1, species are independent from each other. A VR >1 indicates positive correlation, while a VR <1 indicates negative correlation. Statistics W and χ2 test were used to test the significance of the difference between VR value and 1, with a χ2 (0.05, N) < W < χ2 (0.95, N) indicating statistically insignificant di fference between VR value and 1, which suggests the original hypothesis is true. The formulas are listed below [6]:

$$\text{VR} = \frac{\mathbf{S}\_{\text{T}}^{2}}{\delta\_{\text{T}}^{2}} = \frac{\frac{1}{\text{N}} \sum\_{j=1}^{N} \left(\mathbf{T}\_{\text{j}} - \mathbf{t}\right)^{2}}{\sum\_{i=1}^{\text{S}} \left(\mathbf{P}\_{\text{i}} \times \left(1 - \mathbf{P}\_{\text{i}}\right)\right)} \left(\mathbf{P}\_{\text{i}} = \frac{\mathbf{n}\_{\text{i}}}{\mathbf{N}}\right) \tag{1}$$

$$\mathbf{W} = \mathbf{V}\mathbf{R} \times \mathbf{N} \tag{2}$$

where S represents the total species number, N represents the total quadrat number, Tj represents the number of species occurring in the quadrat j, t represents the average number of species occurring in each quadrat, while Pi represents the percentage of the quadrats with the species i occurring in them.

Then, based on 2 × 2 contingency table among species, we used the Yates correlation coe fficient, Ochiai index and principal component analysis (PCA) to reveal the interspecific association between each pair of species. First, Yates correlation coe fficient was processed to test whether the associations of species pairs existed with the predetermined probability. In Yates test, index V > 0 indicates a positive correlation between the two species, while V <0 indicates a negative correlation between species. χ2 is used to show the significance of the correlation. If χ2 < 3.841, the two species are independently distributed, if 3.841 ≤ χ2 < 6.635, significant association is attached, and if χ2 > 6.635, the association between two species is highly significant [6]. Second, Ochiai index was used to measure the degree of the associations. Ochiai index distributed between 0 and 1. A higher OI indicates higher possibility for two species to occur in the same quadrat. The formulas of Yates correlation and Ochiai index are listed below [5,6]:

$$\mathbf{V} = \mathbf{ad} - \mathbf{b}\mathbf{c} \tag{3}$$

$$\chi^2 = \frac{\mathbf{n}(|\mathbf{ad} - \mathbf{bc}| - 0.5\mathbf{n})^2}{(\mathbf{a} + \mathbf{b})(\mathbf{a} + \mathbf{c})(\mathbf{b} + \mathbf{d})(\mathbf{c} + \mathbf{d})} \tag{4}$$

$$\text{OI} = \frac{\text{a}}{\sqrt{(\text{a} + \text{b})(\text{a} + \text{c})}} \tag{5}$$

where n represents the total quadrat number, a represents the number of quadrats with both two species, b and c represent the number of quadrats with only one species, while d represents the number of quadrats with neither of the two species. Finally, PCA analysis was adopted to show the relationship among species and quadrats at the community level.

#### *2.3. Hydraulic Traits Measurement*

Long branches over 50 cm with diameters of approximately 5 to 8 mm were cut by clippers, and all the leaves on the cut branch were retained. The branch samples were placed in black plastic bags, and the cross sections were immediately wrapped in moist tissue in the field to reduce moisture loss. The branch samples were cut twice into 27.4-cm length segments under water in the laboratory [22], removing the excess length from both ends, ensuring shaving o ff more than 5 cm from each end [22,23] and retaining a straight segmen<sup>t</sup> with as few embranchments as possible.

Maximum hydraulic conductivity (Kmax, 10−<sup>4</sup> kg m MPa−<sup>1</sup> s<sup>−</sup>1) was measured for each sample by flushing the stem segments with the 20-mM KCl solution at 0.1 MPa for 30 min and measuring the flow rate of ultrafiltered water passing through the segmen<sup>t</sup> under pressure di fference of around 0.05 MPa [15,24]. All the leaves above the sampled branches were scanned using a portable scanner (Founder MobileO ffice Z6, Beijing, China) and the total leaf area was calculated using MATLAB R2014a (The MathWorks, Natick, MA, USA). Leaf specific conductivity (Kl, 10−<sup>4</sup> kg m<sup>−</sup><sup>1</sup> s<sup>−</sup><sup>1</sup> MPa−1) was then determined as the ratio of Kmax and the total leaf area of the branches. Kl is a synthetic proxy for the balance between the water demand of leaves and the hydraulic transport e fficiency of stems [14]. All the measurements were made within two days after sampling.

The vulnerability of the stems to drought-induced cavitation, i.e., the stem-vulnerability curve, was estimated using the centrifugal force method [25] with a superspeed centrifuge (CTK150R, XiangYi Centrifuge Instrument, Changsha, China). The segmen<sup>t</sup> ends were kept immersed in water during spinning [25]. The Kmax of the stem segments was set as the starting point of the curve. Then, the hydraulic conductance (Kh, 10−<sup>4</sup> kg m MPa−<sup>1</sup> s<sup>−</sup>1) were measured after a series of stepwise increasing xylem tensions caused by di fferent spinning speeds. The vulnerability curve was finally plotted as the variation in the percentage loss of hydraulic conductivity (PLC) with increasing xylem tension. The PLC of the stems was calculated as:

$$\text{PLC} = \left[ (\text{K}\_{\text{max}} - \text{K}\_{\text{h}}) / \text{K}\_{\text{max}} \right] \tag{6}$$

From the vulnerability curves, we can extract some important water potential thresholds, which mark declines in stomatal conductance, hydraulic conductivity and potential for survival. Here, we choose stem water potential at 12% (P12), 50% (P50) and 88% (P88) loss of stem conductivity as the proxies for the early declines in water supply influencing leaf stomatal activities, gas exchange and carbon uptake, the maximum water stress trees can bear naturally and the threshold of irreversible xylem functional damage, respectively [26].

### *2.4. Statistical Analysis*

The stem-vulnerability curves were fitted with a four-parameter Weibull model using SigmaPlot v14.0 (Systat Software, Inc., San Jose, CA, USA). Variance analysis was used to show the di fferences among the hydraulic traits of di fferent species, followed by least significant di fference (LSD) post hoc tests. To reveal the relationship between the interspecific association pattern and the hydraulic-strategy di fferences among species, we calculated the hydraulic di fferences between each pair of species with Kl and P50, including Euclidean, Manhattan, Canberra, Minkowski and maximum distance. If all of the distances show similar relationships with interspecific association indices, we can accept the hypothesis that linkages exist between hydraulic-strategy di fferences and interspecific associations. Correlation tests and linear regression were used further for the hydraulic di fferences and the interspecific association indices, V2012, V2015, OI2012 and OI2015 to map the relationship between them. All the statistical analyses were performed in R software (R Development Core Team, 2009), while figures were drawn with R software and SigmaPlot v14.0.
