Introduction

Drought-tolerant trees and shrubs are critical components of Mediterranean and arid zone ecosystems (Larter et al., 2017; Klein, 2020; Akhmedov et al., 2021). With climate change increasing heatwave and drought intensity and frequency in many of these regions, understanding how vulnerable these species are to environmental stress is important (Willson and Jackson, 2006; Larter et al., 2017; Wagner et al., 2022). Junipers, conifers from the widespread family Cupressaceae, are among the most drought-tolerant genera of trees in the world, and count over 50 species appearing in semiarid regions across the northern hemisphere. Some of these species also grow at remarkably wide seasonal temperature ranges, enduring annual minimum and maximum temperatures of -20 ˚C and +40 ˚C, respectively (Willson and Jackson, 2006; Seim et al., 2016). Considering the ongoing global warming on the one hand, and the increasing occurrence

of extreme climatic events on the other hand (Reyer et al., 2013), junipers can become focal tree species in the near future.

In North America, many juniper species were reported to increase in distribution and abundance in arid and semiarid regions, e.g., across western and central USA (Willson et al., 2008). Among 14 studied juniper species, survival during extreme droughts in southwestern USA and expansion into arid habitats was explained by high xylem resistance to embolism (Willson et al., 2008). For example, Juniperus pinchotii had 50% loss of hydraulic conductance (P₅₀, related to xylem embolism) at a very negative water potential of -14.1 MPa. Among four Arizona Juniper species, divergent resistance to embolism was measured under both freezing and drought backgrounds (Willson and Jackson, 2006).

In drylands from the Mediterranean to central Asia, junipers are widespread, with often disjunct and relictual endemic distributions (Boratynski et al., 2023). As such, they are vulnerable to extinction, especially on the dry margin of their ranges. Drought stress related with climate change was reported to increase in the mountains of central Asia, where junipers growing at lower elevations became sensitive to increasing summer drought (Seim et al., 2016). At higher elevations, juniper growth, previously favored by warm summer temperatures, has become negatively affected by increasing summer aridity. Among Mediterranean junipers, plant ecology and physiology have just recently been investigated (Baquedano and Castillo, 2008; Boratynski et al., 2023; Korakaki et al., 2024). Pollen records from lake sediments in the Levant indicate junipers’ existence along the past millennia in locations where they are missing today (Bar-On et al., 2022).

Recent ecophysiological studies have highlighted that xylem embolism is a key mechanism inducing plant death from individual leaves and branches through to ecosystem scale dieback events (Anderegg et al., 2016; Cardoso et al., 2020; Mantova et al, 2022; McDowell et al., 2022; Sorek et al., 2023). Junipers are known to be an embolism resistant genus, explaining their survival during dieback events that affect co-occurring conifer species (McDowell et al., 2008). However, we still lack a broad understanding of the variation in drought-related functional traits in this genus. P₅₀ is a major index of xylem vulnerability to embolism, indicating the water potential inducing 50% loss of xylem conductivity due to embolism. Variation in P₅₀ across tree species is large, with a P₅₀ being less negative than -4 MPa for many tree species (Choat et al., 2012). Yet, some species stand out with a P₅₀ of -10 MPa or less, notably among the two Cupressaceae subfamilies (Callitroideae in the southern hemisphere and Cupressoideae in the northern hemisphere (Pittermann et al., 2010). Among these, some juniper and callitris species have been identified (Larter, 2016). Interestingly, in terms of global distribution, most of the tree species with high resistance to embolism (P50 of -10 MPa or lower) are from Australia and North America, with none from Eurasia (Willson et al., 2008; Larter et al., 2017).

In conifers, water transport occurs through single-celled tracheids that are interconnected via pit membranes (Hacke and Jansen, 2009; Delzon et al., 2010; Pittermann et al., 2010). Contrary to former belief, embolism resistance is conferred mostly through pit anatomy, rather than tracheid anatomy (Bouche et al., 2014). Pit membranes have two components: the porous margo, allowing water to pass through the membrane, and a central, impermeable, thickening, namely the torus, which isolates gas-filled tracheids. Under drought, tracheids can become air-filled and hence hydraulically dysfunctional when the torus margo fails and an air bubble penetrates a functional conduit, nucleating embolism in the water column. Pit depth and torus thickness correlated positively with embolism resistance in three boreal conifer species (Hacke and Jansen, 2009). Across 15 Cupressaceae species, embolism resistance was most closely correlated with the ratio of the torus to pit aperture diameter, i.e., the torus overlap (Delzon et al., 2010; Pittermann et al., 2010). A large overlap between the torus and the pit aperture allows the torus to tightly seal the pit aperture when needed. This embolism resistance mechanism was corroborated in analyses testing an increasing number of conifer species: among 28 species (Song et al., 2022), among 40 species (Delzon et al., 2010) and among 115 species (Bouche et al., 2014). Within ecotypes of the key Mediterranean species Pinus halepensis, torus overlap and pit aperture size also correlated with P₅₀ variations (David-Schwartz et al., 2016). However, information about Eurasian juniper species is scarce.

Here we studied xylem functional traits of three species of junipers in three distinct field sites, namely Juniperus drupacea (Mt. Hermon, Israel), J. pseudosabina (Zaamin National Park, Uzbekistan), and J. seravschanica (Aman-Kutan mountains, Uzbekistan) with the aim of better understanding their xylem resistance to embolism, and its relationships to xylem anatomy and climatic variables. We then put these results in the context of the global phylogeny of the genus Juniperus. Our hypotheses were that (1) embolism resistance is correlated with aridity, therefore the species growing in drier areas would have the more negative P₅₀, and (2) that embolism resistance would be linked to coordinated changes in xylem anatomy, specifically to pit structure features.

Materials and Methods

Study species, sites, and climate

We studied mature trees of three Eurasian Juniperus species: two from central Asia (J. pseudosabina and J. seravschanica) and one Mediterranean Juniperus species (J. drupacea). The latter is a tree species growing in the eastern Mediterranean region from Greece to Syria and Israel. In this study, we sampled individual trees from a population in Mount Hermon, in northern Israel (Fig. 1; Table 1). J. drupacea trees can reach 20 m in height and are included in the IUCN list of species. J. drupacea study trees were 5-8 m in height and 38-63 cm in diameter at breast height (DBH). J. pseudosabina are shrubs or small trees growing up to 10 m tall at elevations of 2000 to 4000 m in Afghanistan and to the north-east to China and Mongolia. J. seravschanica (J. excelsa subsp. polycarpos) grows up to 20 m tall and occurs in mountain ranges at 1200 m to 3000 m from eastern Turkey to the Caspian Sea, and east into Afghanistan, Pakistan and India.

Table 1: A summary of juniper species and study sites, with major climatic characteristics for each site.

Figure 1: Sampling sites and study trees of the three Eurasian Juniper species. Sites and species match horizontally. For geographical locations see Table 1 and Fig. 2.

J. seravschanica (J. excelsa subsp. polycarpos) grows up to 20 m tall and occurs in mountain ranges at 1200 m to 3000 m from eastern Turkey to the Caspian Sea, and east into Afghanistan, Pakistan and India. The two Asian species were respectively sampled in Zaamin National Park and the Aman-Kutan mountains in Uzbekistan (Fig. 1; Table 1). J. pseudosabina study trees were 5-12 m in height and 70-110 cm in DBH. J. drupacea study trees were 5-12 m in height and 60-110 cm in DBH. For each species, we downloaded occurrences from GBIF (Holstein 2001; https://www.gbif.org), and manually curated them to remove points that landed outside of the natural ranges of the species (Farjon, 2010; Earle, 2024 (conifers.org)). For the remaining points (J. drupacea, n= 75; J. pseudosabina, n=423; J. seravschanica, n=156), we obtained historical climatic data (1970 to 2000) from the worldclim v2.1 layers at 2.5-minute resolution (Fick and Hijmans 2017) using R (R Development Core Team, 2022), which we summarized to obtain median and quantiles of the climate range of each species. Considering drought-induced embolism, we chose to focus on mean annual temperature, mean annual precipitation, and precipitation of the driest quarter. We also calculated the aridity index of the driest quarter, as the ratio between precipitation and potential evapotranspiration of the driest quarter.

Branch sampling

For each species at each site, we selected 5 individual trees, from which we collected several branches in April 2022. Across sites, April marks the end of the wet season, when trees are well-hydrated and frost-free, i.e., when xylem embolism is highly unlikely. We selected straight branch segments from lower, but exposed, parts of the crowns, with no side branches, and with

diameters of around 1 cm and a length of approximately 40 cm. All leaves were removed and the samples were wrapped in damp paper towels, triple-bagged and sent to the Caviplace lab (INRAE, University of Bordeaux, France) for measurement of hydraulic traits.

Hydraulic functional traits

Embolism resistance was characterized by P₅₀, which was calculated based on vulnerability curves measured using the Cochard Cavitron method (Cochard et al., 2013). This technique uses centrifugal force to generate xylem tension in a branch sample and allows the rapid estimation of flow of water through the sample. By measuring xylem conductance at decreasing xylem pressure steps, we can determine the samples’ vulnerability curve, and estimate the xylem pressure inducing 50% loss of hydraulic conductance, P₅₀. K_s was calculated as the hydraulic conductivity measured in the Cavitron at the lowest rotation speed (i.e., the highest water potential, before any conductance is lost to embolism). For each species, we measured 3 branches from each of 5 trees. Since xylem conduits in junipers are tracheids, there is no risk of an “open-vessel” artefact and so-called “r-shaped” curves reported in some angiosperm species (Torres-ruiz et al., 2017). A sigmoid Pammenter model was fit to each individual curves using SAS software (Pammenter and Van der Willigen, 1998). Vulnerability curves, corresponding to percentage loss of hydraulic conductivity (PLC, %) as a function of xylem pressure (P) were fitted using the NLIN procedure in SAS 9.4 (SAS, Cary, NC, USA) based on the following equation:

Eq. (1) $\mathrm{P}\mathrm{L}\mathrm{C}=$100/ 1 + exp (S/25 ∗ (P − P₅₀))

where P₅₀ (MPa) is the xylem pressure/water potential inducing 50% loss of hydraulic conductivity and S (% MPa^-1) is the slope of the vulnerability curve at the inflection point. To place the data in a broader context, we obtained published data for other juniper species (Hammond et al., 2021). All anatomical observations were carried out on material that had previously been used for measuring embolism vulnerability. We selected four individuals for each species, choosing the samples that were closest to the average P₅₀ value. We observed the samples with a benchtop scanning electron microscope at the PHENOBOIS platform of University of Bordeaux_INRAE (PhenomG2 pro; FEI, The Netherlands). Samples 5–8 mm long were cut with a razor blade in a radial direction. After drying for 24 h in an oven at 70 °C, the samples were fixed on stubs, coated with gold using a sputter coater (108 Auto; Cressington, UK) for 40 s at 20 mA, and observed under 5 kV. All anatomical features were based on earlywood, which is responsible for most of the hydraulic conductance. Average values for pit anatomical traits were determined based on at least 20 measurements per sample, made using ImageJ 1.52r (Schindelin et al., 2012). We measured 20 pits per individual, four individuals per species.

In situ xylem water potential

To learn more about the more resistant J. drupacea, further investigation was performed. In September 2022, at the end of the dry season in Mt. Hermon (Israel), a field campaign was carried out. Tips of branches of J. drupacea were wrapped in sealed, tin-foiled plastic bags for 30 min. for equilibration prior to sampling and were then cut and their water potential (WP) was measured in situ. Since transpiration was arrested on these tips, the measured WP could indicate the xylem, and not leaf, WP. The seasonal timing of measurement was chosen to represent the driest conditions on site, with the annual minimum soil moisture (~13 % v/v; Bar-On et al., 2022) and a high vapor pressure deficit (VPD) of >4 kPa. Xylem WP was measured using a pressure chamber (PMS Instruments, Albany, OR, USA) on branch tips from five trees at three cycles, in three cycles: 12:10-12:50; 12:50-13:25; and 13:25-14:00. Xylem WP was measured exclusively for J. drupacea, and not for the other two species.

Data analysis

We tested the effect of species on hydraulic functional traits (P₅₀, torus overlap, and Ks). We applied a generalized linear model approach and compared means using the stepwise multiple comparisons procedure using a Student-Newman-Keuls Test. In addition, linear regressions were calculated for the correlations between P₅₀ and climatic data (temperature and precipitation). To represent the evolutionary relationships of the three study species within the Juniperus genus, we used a published molecular phylogeny (Forest et al., 2018). All statistical tests were performed in R software and in SAS software (SAS® OnDemand for Academics, Cary, NC, USA).

Results

Distribution ranges of the studied species diverged from a relatively small region for Juniperus drupacea (Eastern Mediterranean) to a larger region for J. seravschanica (central and western Asia) and a very large region for J. pseudosabina (central Asia and Northeast into Siberia). Expectedly, range-wide climatic conditions for the two Asian species were larger than for the Mediterranean species, also reflecting the higher climatic variability in continental vs. coastal habitats (Fig. 2). Overall, J. drupacea was confined to a relatively narrow range of mean annual temperatures, which mostly represented the warmer edge of the Asian species. The latter species’ distributions spanned over locations with up to 30 ˚C mean annual temperature difference and mostly overlapped between the species, with J. seravschanica occupying warmer locations and J. pseudosabina occupying colder locations (Fig. 2).

Figure 2: Distributions and climatic envelopes of three Eurasian Juniper species. Distribution ranges were constructed according to occurrence data from GBIF.org (circles). Sampling locations are denoted over the map and within the climatic space of precipitation of the driest quarter and mean annual temperature (inset).

Comparing among the species’ locations by precipitation of the driest quarter (PDQ) showed that most populations inhabited locations with a distinctive dry season (0-50 mm during the driest quarter), while some J. seravschanica and J. pseudosabina populations at higher precipitation locations (50-100 mm during the driest quarter). The specific populations studied here (Table 1) did not necessarily represent the average climatic conditions for each species. For J. drupacea, Mt. Hermon with a mean annual temperature (MAT) of 15 ˚C was on the warmer side of its MAT range from 7 ˚C to 17 ˚C. For J. seravschanica, Aman-Kutan (MAT = 10 ˚C) was on the warmer side of its MAT range from -10 ˚C to 24 ˚C. However, for J. pseudosabina, Zaamin (MAT = 4 ˚C) was in the middle of its MAT range from -14 ˚C to 14 ˚C. Similarly, Mt. Hermon and Aman-Kutan were on the drier side of their species’ PDQ range, while Zaamin was roughly in the middle. Overall, the studied populations spanned between the dry and warm J. drupacea and the cold and wetter J. pseudosabina, with the J. seravschanica population in the middle.

Xylem embolism resistance varied widely across the three studied junipers (Fig. 3). J. drupacea was the most embolism resistant with a P₅₀ = -11.8±0.3 MPa, followed by J. seravschanica (P₅₀ = -8.7±0.8 MPa) and J. pseudosabina (P₅₀ = -6.3±0.6 MPa). Differences in P₅₀ among the three species were significant (P <0.0001). To further investigate the more resistant J. drupacea, xylem water potential (WP) was measured on the same five trees sampled in the field. Overall, WP in J. drupacea at midday at the end of the dry season was -4.85±0.19 MPa.

Among xylem anatomy traits, torus overlap ranged between 0.20 and 0.42 µm across the species, with large intra-species variations. However, there was no difference among species (P = 0.23-1.00 in pairwise comparisons using Kruskal-Wallis test). Therefore, there was no correlation between torus overlap and P₅₀ (Fig. 4).

Maximum xylem specific hydraulic conductivity was also less variable among species than P₅₀, with K_s values ranging from 2.26 10-4 kg m^-1 MPa^-1 s^-1 in J. seravschanica (polycarpos) to 3.62 10-4 4 kg m^-1 MPa^-1 s^-1 in J. pseudosabina. The most embolism resistant species, J. drupacea was intermediate with K_s = 2.86 10-4 kg m^-1 MPa^-1 s^-1. We further plotted P₅₀, torus overlap, and K_s against mean annual temperature and precipitation of the driest quarter for each species and site (Fig. 5). We found that hydraulic trait variation did not relate to the climate of the species, nor to the sampling sites, for both torus overlap and K_s. However, P₅₀aligned well with mean annual temperature at both species and population (site) levels, whereby P₅₀was more negative by ~1 MPa for each 2 ˚C warming (P₅₀= -0.55×MAT -4.98; r² = 0.997; linear regression based on Fig. 5). P₅₀ aligned well also with precipitation of the driest quarter at the site level, whereby P₅₀was more negative by ~1 MPa for each 6 mm drying (P₅₀= 0.17×MAT -11.21; r²= 0.917; linear regression based on Fig. 5). Interestingly, this relationship did not express at the species level, i.e., overall, the most drought tolerant species, J. drupacea, did not inhabit drier sites than the Asian species. Next, P₅₀ was plotted against the aridity index of the driest quarter (Fig. 6). Xylem embolism resistance was higher (P₅₀ was more negative) for species with higher aridity (lower value of aridity index). This strong relationship was evident at both species (P₅₀ = 18.2×AI -12.0; r² = 0.811) and site levels (P₅₀ = 40.57×AI -10.9; r² = 0.835).

Figure 5: Relationships between three hydraulic traits (y-axes) and two climatic variables (x-axes) across three Eurasian Juniper species. Circles and error bars represent averages ± SE for each species based on its entire distribution range (Fig. 2), with other shapes denoting the climatic conditions at the specific study sites.

Finally, we added the newly measured P₅₀ values for the three juniper species (Fig. 3) to P₅₀ values measured previously for 16 juniper species (out of 53 known juniper species worldwide). The three Eurasian juniper species studied here were within the genus P₅₀ range between -6 MPa (J. squamata) and -14 MPa (J. pinchotii). Specifically, only three species had more negative P₅₀ than J. drupacea. Unlike the Mediterranean J. drupacea, those were new world species. There was no evident relationship between the species-specific P₅₀ and aridity at the species’ range. For example, both J. pinchotii and J. drupacea were from less arid sites than those of J. seravschanica (polycarpos) and J. pseudosabina. Notably, the two Asian junipers studied here represented the highest aridity for the genus.

Discussion

In this study we measured hydraulic vulnerability to xylem embolism in three, yet unexplored, Asian juniper species: Juniperus drupacea growing in Israel, and J. pseudosabina and J. seravschanica growing in Uzbekistan. The P₅₀ of the Eastern Mediterranean J. drupacea was -11.8±0.3 MPa, i.e., the most drought-tolerant tree species in Asia, to the best of our knowledge. Although the Asian species J. pseudosabina and J. seravschanica were significantly less drought-resistant, they are still very well adapted to extreme droughts.

Our first hypothesis was that embolism resistance is correlated with aridity, therefore the species growing in drier areas would have the more negative P₅₀. Our results counter this hypothesis, since the Asian species J. pseudosabina and J. seravschanica were less resistant than the Mediterranean J. drupacea (Fig. 3), despite the Asian species growing in very arid environments. Contrary to another Cupressaceae genus with drought affinities, Callitris (Larter et al., 2017), we found little evidence of a strong correlation between embolism resistance and aridity across 19 species of junipers (Fig. 7). How can we reconcile the weak relationships between aridity and embolism resistance in the three Juniperus species? We suggest three potential explanations to this observation. First, climate at the specific Israeli site was warmer and with the driest summer (Fig. 5). In other words, it is not the annual precipitation amount that matters for embolism resistance, but rather the length and intensity of the dry season (Table 1). Indeed, when P₅₀was plotted exclusively against the aridity index of the driest quarter, the expected relationship was observed (Fig. 6). Mt. Hermon in Israel has a long rainless season following heavy snowfalls in winter; both affecting its observed exceptionally low treeline (Bar-On et al., 2023). This explanation agrees with the first in situ measurements of embolism dynamics in trees, performed in Pinus halepensis in a semi-arid forest (Wagner et al., 2022). These measurements showed that embolism occurs during heatwaves, due to the combination of hot and dry conditions (Wagner et al., 2022). A second explanation relates to the exposure of the two Asian species to very cold conditions (Fig. 2). Across four Arizona juniper species, hydraulic vulnerability to a drought-freezing combination was higher than to drought-induced embolism alone (Willson and Jackson, 2006). It is hence possible that freezing (Fig. 2) had a role in decreasing embolism resistance for the Asian species, compared to the Mediterranean species, which is rarely exposed to temperatures <-5 ̊C in Mt. Hermon (Bar-On et al., 2023). The third explanation considers that the Israeli population represents in fact the warm and dry range edge for J. drupacea, far from the species’ average conditions. It is well-known that marginal populations can significantly differ from core populations of a given tree species (Lerner et al., 2023). Specifically for J. drupacea, the Greek populations differ genetically and morphologically from the Asian populations (Boratynski et al., 2023). It is hence possible that P₅₀ of the Greek J. drupacea is less negative than that of J. pseudosabina and J. seravschanica, as expected. However, hydraulic safety (P₅₀) and efficiency were conserved across J. communis ecotypes spanning from southern France to northern Sweden (Unterholzer et al., 2020). Therefore, this explanation is not supported.

Figure 7: Xylem vulnerability to embolism (P₅₀, -MPa) for measured species across the phylogeny of the genus Juniperus. Bar colors denote the 5th percentile of the aridity index of the species distribution (see legend on bottom right corner; high aridity is denoted by low values and red color). The three Eurasian Juniperus species studied here are marked within their respective phylogenetic clades, which are highlighted in color.

To what extent are the studied juniper species under climatic risk? Here, in situ minimum xylem WP of J. drupacea was ~-5.0 MPa. Comparing this value with the species’ vulnerability curve (Fig. 3) indicates zero embolism for these trees during the driest time of the year. This WP value is lower than that of another Mediterranean juniper species, J. phoenicea, which measured -3.5 MPa in dry season midday (Baquedano and Castillo, 2008). WP in J. virginiana was measured down to -6.2 MPa when exposed to a rapid soil dry-down in a greenhouse experiment (Long et al., 2023). In addition, J. drupacea can enforce stomatal closure during the driest quarter, as observed in sap flow behavior in a Greek population (Korakaki et al., 2024). In other words, its hydraulic safety margin is potentially wide, as reported in J. virginiana from the most arid, with P₅₀ as much as 3.25 MPa less negative than minimum shoot water potentials measured (Long et al., 2023). Therefore, the studied J. drupacea trees are probably rather safe at current conditions; suggesting a potential to withstand future drying; highlighting the need to protect these species from direct anthropogenic effects like deforestation (Guo et al., 2022), having higher chance of surviving compared to other tree species. Still, extreme events like severe or long heatwaves with high VPD might push J. drupacea to its hydraulic limit faster than we think (Wagner et al., 2022). Notably, a WP of ~-5.0 MPa as measured here would mean 10-20% loss of conductivity for the Asian species (Fig. 3). Indeed, the J. seravschanica trees studied here in Zaamin National Park were part of a region-wide study reporting increased drought stress observed in junipers of central Asia (Seim et al., 2016).

Our second hypothesis was that embolism resistance would be linked to coordinated changes in xylem anatomy, and specifically to pit structure features. A large overlap between the torus and the pit aperture allows the torus to tightly seal the pit aperture when needed, as shown in numerous conifer species (Delzon et al., 2010, Pittermann et al., 2010, Bouche et al., 2014). Torus overlap values measured here were in the range of 0.20-0.42 µm (Fig. 4) which is typical to many conifers, although slightly on the lower side (Bouche et al., 2014). Within each species, there was rather high variability in P₅₀ (Fig. 3), and moreover, in torus overlap (Fig. 4). Overall, torus overlap alone could not explain the large differences among the juniper species studied here. It is possible, however, that torus porosity, not measured here, played a role, e.g., if torus pores are more prevalent in the xylem of the less embolism-resistant juniper species studied here. In 13 out of 33 conifer species, pores were found, using electron microscopy, within the torus (Jansen et al., 2012). The pores have a plasmodesmatal origin, an average diameter of 51 nm, and arrange in groups. While junipers were not included in that study, three other Cupressaceae species did have pores in the torus (Jansen et al., 2012). Alternatively, could margo traits explain the interspecific P₅₀ differences observed here? Past studies of conifer hydraulics showed mixed results. Margo porosity did not vary systematically with embolism resistance among 15 Cupressaceae species (Pittermann et al., 2010). Among 40 conifer species, highly embolism resistant species exhibited high flexibility of the margo (Delzon et al., 2010). But in a more extensive study of 115 conifer species, the flexibility of the margo did not seem to play a role (Bouche et al., 2014). Last, we also tested the relationship between xylem hydraulic conductivity and climatic variables, which were weak (Fig. 5). Studying 28 conifer species (two juniper species included) grown in a 50-year-old common garden experiment, hydraulic conductivity was only negatively correlated with tracheid cell wall thickness (Song et al., 2022). It remains to be tested if this tracheid trait played a role here.

This study was motivated by the remarkable climatic adaptability of Eurasian junipers on the one hand, and the pressing need to find tree species that would stand future climate extremes on the other hand. Among conifer species with P₅₀ more negative than -10 MPa, only two are outside the Cupressaceae; a single family that includes 19 such embolism-resistant species (Larter, 2016). Among these 19 Cupressaceae species, nine belong to the Australian genus Callitris, where P₅₀ variation was driven by aridity up to its most extreme value of -18.8 MPa, in Callitris tuberculata (Larter et al., 2015; Larter et al., 2017). Here we measured P₅₀ of -11.8±0.3 MPa in J. drupacea growing in Israel, indicating a good potential for survival under climate change. Junipers, including the J. pseudosabina and J. seravschanica studied here, dominate the landscape under the continental climate of central Asia. In North America, juniper species expand into arid habitats across the western and central USA (Willson et al., 2008). What are the chances for such a scenario to happen in Eurasia? For example, Pinus halepensis populations were projected to gradually move northward with climate change (Patsiou et al., 2020). The ability of the three juniper species to follow their climatic niche will likely be limited by elevation and geographical barriers. For example, in Israel, junipers only grow in two mountain sites, while lower regions are dominated by pines and oaks. But pines and oaks are presumably less drought-resistant than junipers, and hence, if we use junipers in current afforestation programs, we might be able to increase forest sustainability into the 21st century.

Acknowledgements

Adriano Sofo (University of Basilicata, Italy), Assaf Yaakobi, Yedidya Ben-Eliyahu, and Inon Aizik (The Weizmann Tree Lab) are acknowledged for their help in the Mt. Hermon field campaigns. The authors thank Yael Wagner (The Weizmann Tree Lab) who commented on an earlier version of this manuscript. Marylou Mantova (University of Florida), Emilie Isasa (University of Toronto) and an additional anonymous reviewer are acknowledged for their valuable comments during peer review. TK thanks the Edith and Nathan Goldenberg Career Development Chair; Mary and Tom Beck-Canadian Center for Alternative Energy Research; Larson Charitable Foundation New Scientist Fund; Yotam Project; Dana and Yossie Hollander; the Estate of Emile Mimran; and the Estate of Helen Nichunsky.

Author contributions

The study was initiated by TK, and most analyses were done by ML. Field samples were collected by AA and TK. Cavitron measurements were performed by CP and SD. Most figures were prepared by ML. The manuscript was written by TK, with contributions from all authors.