Incorporating Ecosystem Functional Diversity into Geographic Conservation Priorities Using Remotely Sensed Ecosystem Functional Types

Conservation biology must set geographic conservation priorities not only based on the compositional or structural but also on the functional dimensions of biodiversity. However, assessing functional diversity is challenging at the regional scale. We propose the use of satellite-derived Ecosystem Functional Types (EFTs), defined here as patches of land surface that share similar primary production dynamics, to incorporate such aspects of ecosystem functional diversity into the selection of protected areas. We applied the EFT approach to the Baja California Peninsula, Mexico, to characterize the regional heterogeneity of primary production dynamics in terms of EFTs; to set conservation priorities based on EFT richness and rarity; and to explore whether such EFT-based conservation priorities were consistent with and/or complementary to previous assessments focused on biodiversity composition and structure. EFTs were identified based on three ecosystem functional attributes derived from seasonal dynamics of the Enhanced Vegetation Index: the annual mean (proxy of primary production), the seasonal coefficient of variation (descriptor of seasonality), and the date of maximum (indicator of phenology). EFT-based priorities identified 26% of the peninsula as being of extreme or high priority and reinforced the value of the ecosystem functional diversity of areas already prioritized by traditional conservation assessments. In addition, our study revealed that biodiversity composition- and structure-based assessments had not identified the full range of important areas for EFT diversity and tended to better capture areas of high EFT rarity than those of high EFT richness. Our EFT-based assessment demonstrates how remotely sensed regional heterogeneity in ecosystem functions could reinforce and complement traditional conservation priority setting.

conservation priorities based on EFT richness and rarity; and to explore whether such EFT-based conservation priorities were consistent with and/or complementary to previous assessments focused on biodiversity composition and structure. EFTs were identified based on three ecosystem functional attributes derived from seasonal dynamics of the Enhanced Vegetation Index: the annual mean (proxy of primary production), the seasonal coefficient of variation (descriptor of seasonality), and the date of maximum (indicator of phenology). EFT-based priorities identified 26% of the peninsula as being of extreme or high priority and reinforced the value of the ecosystem functional diversity of areas already prioritized by traditional conservation assessments. In addition, our study revealed that biodiversity composition-and structure-based assessments had not identified the full range of important areas for EFT diversity and tended to better capture areas of high EFT rarity than those of high EFT richness. Our EFT-based assessment demonstrates how remotely sensed regional heterogeneity in ecosystem functions could reinforce and complement traditional conservation priority setting.

INTRODUCTION
Contemporary conservation planning faces the challenge of safeguarding the ecological processes required for the persistence of biodiversity over time (GBO4 2014) and for the supply of ecosystem services to people (Costanza 2012). To this end, protected areas must represent the most important areas for in situ global conservation effort (Watson and others 2014). Initially, opportunism and aesthetic values drove protected area creation (Palomo and others 2014;Baldi and others 2017). More recently, reserve selection under systematic conservation approaches (Margules and Pressey 2000) has mainly relied on compositional and structural dimensions of biodiversity (for example, Rodrigues and others 2004;Lamoreux and others 2006). However, despite important advances to the design of more comprehensive protected area networks, geographic conservation priorities have seldom considered heterogeneity in ecosystem functions (Callicott and others 1999;Mace 2014;Turner and Gardner 2015). The need for more representative global protected area networks (Visconti and others 2019) that account for the three dimensions of biodiversity (composition, structure, and function ;Noss 1990) could greatly benefit from the explicit inclusion of ecosystem functions and processes that support biodiversity and ecosystem services (Meyer 1997;Cabello and others 2012;Pettorelli and others 2018;Lecina-Díaz and others 2019).
Functional diversity, ranging from gene expression to landscape processes, is an important biodiversity component to be assessed by conservation programs, as it links biological diversity with ecosystem functioning (Cadotte and others 2011;Díaz and others 2007;Chapin and others 2010; Asner and others 2017), services (Balvanera and others 2006;Duncan and others 2015), and resilience (Mouchet and others 2010). Functional diversity estimates have been made by grouping species into functional types based on structural (for example, shrubs, trees, and so on), phylogenetic (for example, Coniferae, Poaceae, and so on), or metabolic strategies (for example, C3, C4, and so on) related to meaningful biological processes Garnier 2002, Lavorel andothers 2007) or by using morphofunctional species traits (Malaterre and others 2019). However, the capacity for species functional types and traits to represent variations in ecosystem functional properties at regional scales remains a challenge (Wright and others 2006;Pasari and  Understanding the causes and consequences of spatial heterogeneity in ecosystem functions could help protect the species and communities that they support (Meyer 1997;Lovett and others 2005;Turner and Gardner 2015) and elucidate the links between ecosystem multifunctionality, ecosystem services (Manning and others 2018) and ecological stability (Oliver and others 2015). Environmental heterogeneity is a universal driver of taxonomic, phylogenetic, and functional diversity (Stein and others 2014;Stark and others 2017;Bergholz and others 2017). However, although conserving biophysical setting variability has been suggested to preserve biodiversity against rapid environmental change (for example, Lawler and others 2015; Littlefield and others 2019), variation in ecosystem functions has received less attention (Lovett and others 2005). Developing feasible approaches to understand and account for heterogeneity in ecosystem functions could complement traditional priority settings to achieve the holistic goal of protecting all biodiversity facets.
Satellite remote sensing can guide conservation actions based on the characterization of functional diversity not only at the species trait level (Jetz and others 2016) but also at the ecosystem level (Cabello and others 2012; Alcaraz-Segura and others 2013; Asner and others 2017; Gamon and others 2019). First, satellite-derived descriptors of ecosystem functions can be relevant as essential biodiversity variables (EBVs, others 2016, 2018; Alcaraz-Segura and others 2017). For example, spectral indices are linked to key ecosystem functional descriptors such as primary production, evapotranspiration, surface temperature, and albedo (Paruelo and others 1997; Ferná ndez and others 2010; Lee and others 2013) (Ta-Incorporating Ecosystem Functional Diversity in Conservation Priorities ble 1-steps 1 and 2). Second, with these descriptors, it is possible to identify and map areas sharing similar dynamics of matter and energy exchange between biota and physical environments based on so-called satellite-derived Ecosystem Functional Types (EFTs) (Paruelo and others 2001;others 2006, 2013).
As highlighted by Mucina (2019), EFTs could represent ''the first serious attempt to group ecosystems (at large scales) on the basis of shared functional behavior.'' EFTs group ecosystems on the basis of shared ecosystem functions without prior knowledge of vegetation types or canopy structure (Ivits and others 2013). As species can be grouped into plant functional types based on common morphofunctional traits to derive ecological properties at higher biological levels (that is, a bottom-up strategy), ecosystems can be grouped into EFTs to directly map processes and functions at larger scales (that is, a top-down approach) (Alcaraz-Segura and others 2006). EFTs follow a holistic approach (Naeem 1998(Naeem , 2002Loreau 2008) to measure the overall performance of an ecosystem (see the review in Jax 2010). EFTs capture heterogeneity in ecosystem functions (for example, primary production, evapotranspiration, or disturbance dynamics) and provide complementary information to other metrics such as those of vegetation structure and species composition to improve our understanding of the multidimensional nature of biodiversity (Noss 1990). EFTs have already been used to characterize the spatial heterogeneity of ecosystem functioning at the global (Ivits and others 2013), regional (Paruelo and others 2001;Alcaraz-Segura and others 2006;Lara and others 2017), and protected area scales (Ferná ndez and others 2010; Cabello and others 2013).
In this study, we propose the use of Ecosystem Functional Types (EFTs), defined here as patches of land surface that share similar primary production dynamics (that is, productivity, seasonality, and phenology, Figure 1), to incorporate the spatiotemporal heterogeneity of a focal ecosystem function into geographic conservation priorities (conceptual workflow shown in Table 1). As a proof of concept, we applied the EFT approach to the Baja California Peninsula (Mexico): (1) to characterize the regional heterogeneity of primary production dynamics using EFTs; (2) to prioritize areas for conservation based on their EFT diversity (EFT richness and rarity); and (3) to explore whether such EFT-based priorities were congruent with and/or complementary to previous expert and systematic conservation-based assessments mainly focused on biodiversity composition and structure.

Study Area
We chose the Baja California Peninsula as study area ( Figure 2A) because it has high environmental heterogeneity, low human influence, a large proportion of protected land (40%) (Table S1) and because two geographic priority assessments have been conducted on the area mainly based on biodiversity composition and structure (Arriaga and others 2000;Koleff and others 2009). The peninsula covers a Mediterranean desert tropical climatic transition area positioned along a 1400 km latitudinal gradient from 35ºN to 23ºN (Gonzá lez-Abraham and others 2010). The Mediterranean Region (NW) is characterized by annual mean temperatures between 8 and 21°C, dry summers and mild wet winters with annual rainfall levels ranging from 100 to 200 mm at sea level to 500-700 mm in the highest mountains (3100 m) (Peinado and others 2011). The Desert Region extends from NE to S and is characterized by temperatures ranging from 20 to 25°C, and very low annual rainfall (100-200 mm) concentrated in sporadic events that shift from the winter in the N to the summer in the S (Hastings and Turner 1965). The Tropical Region at the southern tip is warm year round (15-24°C) and characterized by a nine-month dry season (November-July) followed by the tropical cyclone and storm rains with annual rainfall levels ranging from 200 mm at sea level to 700 mm in the highest mountains (at 2090 m) (Peinado and others 2011).

Identifying Ecosystem Functional Types
Regional heterogeneity in ecosystem functions was characterized by identifying Ecosystem Functional Types (EFTs) based on the seasonal dynamics of carbon gains following Alcaraz-Segura and others (2013). We focused on primary production because it is an integrative component of ecosystem functioning ( Step 1. To identify the targeted functional facets of biodiversity to be considered, for example, ecosystem primary production as an essential biodiversity variable Conservation planning based on functional dimensions of biodiversity is needed (Noss 1990) but scarce (Cabello and others 2012). Some facets of ecosystem functioning are more essential to biodiversity and ecosystem services, offer more available information for inventorying and monitoring, and are more relevant to particular conservation goals than others We chose primary production, as it is the most integrative indicator of ecosystem functioning (Virginia and Wall 2013) Step 2. To choose surrogates for targeted functional facets, for example, remotely sensed vegetation indices Direct measurements of biodiversity variables are usually costly. Satellite images of the Earth can be considered biological datasets (Geller and others 2017). Image pixels are sampling plots whose spectral information offers indirect, cost-effective estimates of matter and energy exchanges between the land surface and the atmosphere, which support ecosystem functions and services We used the Enhanced Vegetation Index (EVI) to estimate photosynthetically active radiation absorbed by vegetation (based on the Monteith Model, 1972) Step 3. To identify simple and biologically meaningful metrics of the ecosystem functioning surrogates, for example, descriptors of the amount and timing of carbon gain dynamics The dynamics of ecosystem functioning are tracked through full time-series of essential variables. Synthesizing and capturing most of the variance of these time-series into a few easy to interpret metrics reduce complexity, ease interpretability, and promote the metrics standard use We identified three metrics capturing most of the variance in the EVI seasonal dynamics (Ecosystem Functional Attributes, EFAs): annual production, seasonality, and phenology. We parameterized yearly seasonal dynamics of the EVI for three EFAs: the annual EVI mean, seasonal EVI coefficient of variation, and the date of the maximum EVI Step 4. To group patches of the land surface with similar functional behaviors by classifying continuous metrics into discrete units, for example, Ecosystem Functional Types (EFTs) Functional classifications synthesize continuous large-scale ecological gradients into discrete mapping units in relation to common ecosystem functions and processes. Discrete mapping units characterize ecosystem diversity at the regional scale and are needed for management and decision-making such as in systematic conservation planning To integrate patterns of productivity, seasonality, and phenology into a single map, we divided the range of values of each EFA into four intervals (quartiles), creating a potential number of 64 EFTs (4 9 4 9 4) Step 5. To select criteria for assessing ecosystem functional diversity at the regional scale, for example, EFT richness and rarity Measurements of all biodiversity facets are not possible given the complex, multidimensional, and hierarchical nature of biodiversity (Noss 1990). Biodiversity indices such as richness and rarity are easy to interpret, relevant, and objective criteria frequently used in conservation assessments We calculated EFT richness by counting the number of EFTs in a slicing window. EFT rarity was calculated as the relative extension of each EFT compared to the most abundant EFT primary production), the EVI seasonal coefficient of variation (EVI sCV; a descriptor of seasonality), and the date of the maximum EVI (EVI DMAX; an indicator of phenology). The three metrics capture most of the variance in EVI seasonal dynamics into three meaningful metrics that facilitate ecological . The X axis corresponds to months, and the Y axis corresponds to EVI values. EFAs include: the annual EVI mean, an estimator of annual productivity (EVI mean); the seasonal EVI coefficient of variation (EVI sCV), that is, differences between minimum and maximum EVI values, as a descriptor of seasonality; and the date of the maximum EVI (EVI DMAX) as a phenological indicator of the growing season.

Table 1. continued
What is the goal of this step? Why is it needed? How did we complete it?
Step 6. To set geographic conservation priorities that capture areas of high ecosystem functional diversity, for example, areas of high EFT richness and rarity Landscapes of high heterogeneity in ecosystem functions are prone to contain multiple ecosystem metabolic and evolutionary pathways. Multifunctional landscapes provide more diverse ecosystem services (Manning and others 2018), and functional diversity confers ecological stability (resistance and resilience) We identified areas of the highest (extreme and high) conservation priority as those ones with high EFT richness and high EFT composition rarity Step 7. To compare priorities based on ecosystem functional diversity with independent assessments, for example, complementarity and consistency between EFT-based priorities and previous assessments focused on composition and structure Priorities based on ecosystem functioning can converge with independent priorities focused on biodiversity composition and structure so that they reinforce each other. Priorities can also be complementary, supporting decision-making by offering supplementary arguments for the holistic conservation of biodiversity We integrated the three approaches into two synthetic maps: consistency and complementarity. To visualize agreement and disagreement between and among approaches, we used Venn diagrams interpretation (Paruelo and others 2001;Alcaraz-Segura and others 2006). To derive EFT classes from EFAs, the range of values of each EFA was divided into four intervals that were then combined, generating a potential number of (4 9 4 9 4) 64 EFTs ( Figure S1D and S2). We used this classification method with fixed boundaries between classes to maximize the biological interpretability of EFTs and to apply the same classification rules to each year. This way, the classification can be used to track interannual changes in spatial heterogeneity of ecosystem functions (Littlefield and others 2019). As for DMAX since we wanted to maintain its ecological sense in our final classification (that is, the timing or phenology of the interception of radiation by vegetation), the four intervals agree with the four seasons of the year: spring (April-June), summer (July-September), autumn (October-December), and winter (January-March). For EVI_Mean and EVI_sCV, we extracted the first, second, and third quartiles (that is, the 25th, 50th, and 75th percentiles, respectively) for each year. Then, we calculated the interannual means of the quartiles (average of the 17-year period), which were used as thresholds among classes ( Figure S1D). The four intervals created for each variable produced a relatively low number of potential classes (64) and maintained the EFAs spatial patterns ( Figure S1 and S2).
To code EFTs, we used two letters and a number ( Figure S1D): the majuscule indicates primary production (EVI mean) increasing from A to D; the minuscule represents seasonality (EVI sCV) decreasing from a to d; and numbers are a phenological indicator of the growing season (EVI DMAX): 1-spring, 2-summer, 3-autumn, and 4-winter. To summarize ecosystem function patterns of the 2001-2017 period, for each pixel we calculated the most common EFT (the mode) from the 17 annual EFT maps (Table 1-step 4). We excluded from analyses pixels with human influence according to the human footprint index (HF > 0.5) (Gonzá lez-Abraham and others 2015) and those including anthropogenic land uses in the 2017-updated land cover map (INEGI 2017).

Mapping Geographic Conservation Priorities from EFT Richness and Rarity
To identify geographic conservation priorities based on spatial heterogeneity in our focal ecosystem function (that is, primary production dynamics), we derived two diversity metrics from the EFT map: EFT richness and EFT rarity (Table 1-step 5). Both richness and rarity are indices that are easy to interpret, objective, and commonly used in ecology and conservation (Perrin and Waldren 2020). Richness measures different types of entities in a sample. EFT richness was calculated by counting the number of different EFTs within an 8 9 8 pixel sliding window across the study area, serving as an indicator of spatial heterogeneity in primary production dynamics. From the EFT richness of each year, we obtained the interannual average of EFT richness (Alcaraz-Segura and others 2013). We chose this window size because it includes 64 pixels, which is the potential maximum number of EFTs in our classification. The use of smaller window sizes resulted in many windows reaching the maximum number of classes while larger windows produced too coarse outputs (Appendix 5).
Rarity has also been a central focus in conservation (Soulé 1986). According to its abundancebased definition, rarity refers to how frequently an entity is found within an area (Kondratyeva and others 2019). The rarity of each EFT was used as an indicator of distinctive characteristics (that is, singularity) in primary production dynamics, which are likely to exhibit unique biodiversity features with conservation interest (Meyer 1997). EFT rarity was calculated as the extension of each EFT relative to the most abundant EFT throughout the peninsula (Eq. 1) (Cabello and others 2013).

Rarity of EFT
where Area_EFT max is the area occupied by the most abundant EFT throughout the study area and Area_EFT i is the area of the i EFT evaluated with i ranging from 1 (Aa1) to 64 (Dd4). An average rarity map for all years was obtained, serving as our estimate of regional patterns of ecosystem functional singularity.
To determine EFT-based geographic conservation priorities, we searched for areas of high EFT richness and rarity (Table 1-step 6). First, we stretched (by spatial averaging) the spatial resolution of the EFT rarity map (230 m/pixel) to match the EFT richness map resolution (that is, an aggregated value for 8 9 8-pixel windows). Second, the range of values of both priority criteria variables was divided into four intervals using quartiles. Third, a decision matrix with 4 9 4 = 16 possible combinations of richness and rarity levels was produced. Finally, the 16 combinations of richness and rarity levels were grouped into four final priority categories ( Figure 3A): extreme, high, moderate, and low for combinations that summed to 8, 7, 6, and 5, Incorporating Ecosystem Functional Diversity in Conservation Priorities respectively. Combinations with lower sums were deemed not a priority.

Assessment of Spatial Congruence and Complementarity Between the Functional Approach and Previous Assessments
We explored the congruence and complementarity between the EFT-based geographic conservation priorities and two previous assessments based on compositional, structural, and threat features of biodiversity (Table 1-step 7). The ''systematic conservation'' study by Koleff and others (2009) used robust spatial analysis algorithms in a grid to identify four levels of ''Priority Sites to Conserve'' based on diversity of and threats to vertebrates, plants, and vegetation types. The ''expert-based'' study by Arriaga and others (2000) identified ''Terrestrial Priority Regions'' through qualitative expert workshops that combined multiple biological criteria (that is, species richness and endemicity, centers of diversification and domestication, vegetation types, and so on) with criteria for threats and opportunity (that is, habitat loss and fragmentation, unsustainable management, threatened species, and so on).
For the congruence analysis, we overlapped the three approaches at an 8 9 8-pixel window resolution into two synthetic maps: one that integrated congruence between the approaches (where priorities agreed) ( Figure 3C) and another that revealed complementarity (where priorities did not agree) ( Figure 3D). Congruence with other approaches was defined as the existence of a spatial overlap between EFT-based priorities and one or both of the other approaches. Complementarity with other approaches was defined as the existence of spatial discordance between EFT-based priorities and the previous priorities.
To visualize agreement and disagreement between approaches, we used Venn diagrams and the Sorensen-Dice similarity index (Figure 4). Additionally, to show how our EFT-based approach provides useful and orthogonal conservation priority information relative to traditional approaches, we explored the characteristics of congruent and complementary areas among approaches in terms of EFT richness and rarity ( Figure 5) and of EFAs and EFT frequency (Appendix 4).

Sensitivity Analyses
To assess the effect of the sliding window size (Appendix 5), we calculated EFT richness, rarity, and priorities for double and triple window side lengths (that is, 8 9 8-, 16 9 16-, and 24 9 24pixels). To assess the effect of the number of EFT classes considered (Appendix 6), we calculated EFT richness, rarity, and priorities by reducing the number of EFT classes by 86% (8 classes) and 58% (27 classes). Both effects were assessed three ways: from Pearson correlations between different output maps, from the spatial consistency among different output maps, and from the total percentage of the peninsula prioritized by each output map. Finally, we also assessed the effects of different thresholds of EFT richness and rarity on congruence and complementarity between approaches by means of the Sorensen-Dice similarity F-1 index (Appendix 7).

Regional Patterns of Focal Ecosystem Function by Means of EFTs
All 64 potential EFTs were identified in the Baja California Peninsula ( Figure 2B) and exhibited contrasting distributions across the three main ecoregions of the peninsula (Figure 2A; Gonzá lez-Abraham and others 2010). In the Mediterranean Region to the northwest, EFTs were characterized by moderate-high primary production, moderatelow seasonality, and spring EVI maxima (Figures S1 and S2). The central and northeastern Desert Region was characterized by EFTs with low primary production, low to moderate seasonality, and winter EVI maxima in the center and in various seasons in the northeast. The southern part of the Desert Region was characterized by slightly higher level of primary production and seasonality and by summer-autumn EVI maxima. The Tropical Region in the south was characterized by high levels of primary production and seasonality and by summer EVI maxima ( Figures S1 and S2).

Conservation Priorities Based on EFT Richness and Rarity
EFT richness and rarity ( Figure 2C, D) varied across the peninsula following a combination of latitudinal, longitudinal, and topographical gradients ( Figure 2A) and were found to be partially correlated. Areas of high EFT rarity ranged from low to high EFT richness while areas of high EFT richness always showed high levels of EFT rarity (Figure S3). EFT richness levels ranged from 1 to 26 EFTs per sliding window. Most windows of the highest EFT richness (12-26 EFTs) occurred north   Figure 2A, B). (D) Complementarity among geographic conservation priorities obtained by the three approaches (disagreement between Figure 2A, B). White areas were pixels where none of the categories on the map were satisfied. of 30°N in the Mediterranean Region, where climatic gradients translate into high heterogeneity in EFAs, especially along the mountain divide (Figure 2A). An intermittent fringe of high EFT richness was also found along mountains from the southern San Felipe Desert to the center of the Desert Region (from 31º N to 27°N) and continued southwards along the western desert piedmonts and around wetlands and mangroves (from 27º N to 24°N). Moderate EFT richness (7-12 EFTs) was observed in the Mediterranean mountains, San Felipe Desert, Colorado Delta, mid-mountains along the Gulf Coast (from 26º N to 30º N), and desert areas of the central peninsula. Extensive areas with the lowest EFT richness (1-3 EFTs) were found in plains and piedmonts of the Central and Vizcaíno Deserts, along the southern desert mountains (Giganta Ranges), and in the Tropical Region.
EFT rarity gradients were more pronounced than EFT richness gradients ( Figure 2D). The highest rarity (0.8-0.9) occurred in the northwestern quarter of the peninsula above 30°N (Mediterranean Region), the central eastern desert transition, and around wetlands and mangroves. The Pacific northwestern Central and Vizcaíno Deserts (north from 27º N) showed low rarity (0.4-0.7). The lowest rarity (below 0.3) occurred along Giganta Ranges and in the Tropical Region (south of 28º N). This region, dominated by drought deciduous plant functional types, was mostly occupied by one extensive EFT with high productivity and seasonality and by summer EVI maxima (Da2).
The highest priority areas were found in heterogeneous areas across the Mediterranean Region, the northern and central eastern Desert Region, and around wetlands and mangroves ( Figure 3A). Extreme priority areas occupied 9.6% of the peninsula surrounded by areas of high (16.4%), moderate (18%), and low priority (16.6%). The rest of the peninsula (39.5%) was classified as a nonpriority area for EFT diversity.

EFT-Based Priorities Versus Composition and Structure-Based Approaches
EFT-based conservation priorities partially aligned with other approaches (Figure 3A, B). Five percent of the peninsula was considered to be of the highest priority for all three approaches (Figure 4) and mainly the Mediterranean Region along mountain tops and the Desert Region in isolated areas of mountains, wetlands, and mangroves ( Figure 3C). An additional 14% of the peninsula was prioritized by the EFT-based approach and by either the systematic conservation approach (7%) or expertbased approach (7%) (Figs. 3C and 4).
The EFT-based approach also revealed complementary areas not prioritized by the two previous approaches (7% of the peninsula; Figure 4). These areas were mainly located along mountainsides and piedmonts with riverine systems in the Desert Region: the San Felipe Desert to the northeast, the Gulf coastal desert in the center of the peninsula, and scattered areas along the southern desert (north and south of Magdalena Bay) ( Figure 3C). Conversely, some areas (5% of the peninsula) were prioritized by the two previous approaches but not by the EFT-based approach. This occurred mainly in the Mediterranean mid-mountains and in coastal plains of the central and southwestern deserts (Figs. 3 and 4).
EFAs and EFTs slightly differed among areas prioritized by each approach (Figures S6 and S7). Expert-based priorities (Arriaga and others 2000) were biased toward EFTs with less primary production than the other approaches. Systematic conservation priorities (Koleff and others 2009) were biased toward EFTs with higher primary production than the other approaches. In contrast, EFT-based priorities showed a more unbiased distribution of EFA values and EFT compositions than previous priorities (Figures S4 and S5).
EFT richness and EFT rarity were found to be much higher within areas consistently prioritized by the three approaches (6 and 26 EFTs per 8 9 8pixel sliding window with most richness values from 8 to 13 EFTs and with EFT rarity ranging from 0.8 to 1) than within areas prioritized by only one of the three approaches ( Figure 5). In contrast, areas prioritized only by traditional approaches were biased toward areas of low EFT richness (less than 6) but maintained moderate to high values of EFT rarity (greater than 0.5), especially in the systematic conservation approach. Indeed, despite systematic conservation planning and the expertbased approach performing very similarly in capturing EFT richness, systematic conservation planning tended to better represent areas of high EFT rarity ( Figure S6).

Robustness Against Window Size, the Number of Classes, and Priority Thresholds
The sensitivity analyses revealed that our approach to setting priorities was robust against changes in window size and the number of EFT classes (Appendixes 5 and 6). Correlations of EFT richness and EFT rarity across the 8 9 8-pixel window and coarser window sizes ranged from 0.84 to 0.98 (Table S2) and those between the 64 EFT classes and fewer classes ranged from 0.67 to 0.94 (Table S3). Regional patterns of EFT richness, rarity, and priority were largely consistent across window sizes (85% agreement among final priority maps, Figures S7 and S8) and the number of EFT classes (70% agreement among final priority maps, Figures S9 and S10). EFT-based priorities always exhibited more similarities with the more robust systematic conservation approach than with the qualitative expert-based approach independent of thresholds of EFT richness and rarity used ( Figure S11).

DISCUSSION
Contemporary conservation paradigms aim to maintain all biodiversity dimensions (Noss 1990), including the ecological processes and functions Figure 5. Congruence and complementarity among the three approaches to capture Ecosystem Functional Type (EFT) diversity. Density histograms show the frequency of EFT richness (A) and rarity (B) in areas consistently prioritized by the three approaches (''congruence across all priorities'') and in areas exclusively prioritized by one of the approaches but not by the others (''complementarity across priorities''). Our EFT-based approach focuses on two aspects of ecosystem functional diversity (EFT richness and rarity, i.e. heterogeneity and singularity) while the two other approaches focus on biodiversity composition, structure, and threats based on expert knowledge (Arriaga and others 2000) and systematic conservation planning (Koleff and others 2009). that sustain ecosystem services (Meyer 1997;Mace 2014;Prober and others 2019). In this study, we used satellite-derived EFTs (Paruelo and others 2001), defined here as functionally homogeneous land patches in terms of primary production dynamics, to describe spatial patterns of a focal ecosystem function. We used this focal ecosystem function because it is considered to be an integrative surrogate of stocks and fluxes of matter and energy derived from biological activity (Virginia and Wall 2013) and can be easily characterized by remote sensing. In essence, EFTs allowed us to map the spatial patterns of two indicators of ecosystem functional diversity at the regional scale, that is, EFT richness and EFT rarity. From these patterns, we set geographic conservation priorities based on an ecosystem function that helped us identify important areas for the three dimensions of biodiversity (structure, composition, and function) and highlight complementary areas for this ecosystem function not prioritized by traditional approaches.

Regional Patterns of Ecosystem Functional Heterogeneity
Maps of EFAs, EFTs, and EFT richness and rarity offer a characterization of ecosystem functional heterogeneity of the Baja California Peninsula. This heterogeneity results from a combination of latitudinal, longitudinal and topographic gradients. Such gradients determine strong differences across the peninsula in terms of seasonal dynamics of radiation, temperature, precipitation, evapotranspiration, and vegetation access to groundwater (Peinado and others 2011;Villarreal and others 2016) and have been identified as important for plant diversity (Garcillá n and Ezcurra 2003) and endemism (Riemann and Exequiel 2007).
The highest levels of EFT richness were found where topography and spatiotemporal climate variability maximize ecosystem functional heterogeneity, mainly along mountains and piedmonts of the Mediterranean and Desert Regions. The Mediterranean climate imposes two limitations on plant growth: summer drought and winter cold temperatures (Hastings and Turner 1965). These limiting factors of plant growth are strongly heterogenized by steep altitudinal and orientation gradients (Peinado and others 2011). In the Desert Region, latitude, orientation, and access to groundwater impose varying constraints on plant growth. Such constraints include the latitudinal change in the proportion of winter and summer rains; the influence of coastal fog (Webb and Starr 2015); and the occurrence of shallow aquifers, gullies and dry arroyos embedded within a dryland matrix (Leó n de la Luz and others 2015). Such high contrasts in ecosystem functions between the regional landscape matrix and its embedded ecosystems (that is, less water-limited EFTs within a matrix of dryland EFTs) enhance ecological processes of the lateral transfer of matter and energy (Turner and Gardner 2015). For these reasons, despite being a desert, such high heterogeneity in environmental factors renders the Desert Region very diverse in EFTs, a pattern also found for plant functional types and plant communities (Webb and Turner 2015).
The lowest levels of EFT richness were found in the tropics due to wetter and highly consistent tropical climatic conditions that homogenize vegetation (Peinado and others 2011). In the Tropical Region, strong precipitation seasonality (summerautumn tropical rains followed by a nine-month drought) concentrates the growing season following the cyclone season (Leó n de la Luz and others 2000). This high level of seasonality neutralizes even the altitudinal heterogeneity of the mountains, resulting in a spatial homogenization of primary production dynamics throughout the region. Such low EFT richness agrees with high similarities in vegetation composition along all topographic gradients, dominated by a few dry deciduous shrubs and trees (Rascó n-Ayala and others 2018). Such an effect penetrates northwards along the Giganta Ranges with similar vegetation types to the Tropical Region (Gonzá lez-Abraham and others 2010). In addition, very low EFT richness extended northwards along Central and Vizcaíno desert plains and piedmonts. EFT richness in these piedmonts, where energy and water were decoupled (winter rains dominate the Pacific northwestern Central and Vizcaíno deserts, north from 27º N), was lower than in piedmonts where energy and water were coupled (summer rains dominate the southern half of the peninsula and San Felipe Desert to the northeast; Figure S1C).
EFT rarity was found to be associated with latitude, altitude, and the presence of contrasting ecological conditions. The highest EFT rarity of the Mediterranean Region and San Felipe Desert was found to be associated with winter precipitation, which creates a rare phenological pattern in the peninsula (Peinado and others 2011) together with the longitudinal gradient and topographical heterogeneity (for example, the only region with areas showing EVI maxima in all seasons). In the ecological transitional zone of the center of the peninsula (28-29º N), the combined influence of summer tropical storms from the south and autumn-to-spring fronts from the north (Gonzá lez-

Incorporating Ecosystem Functional Diversity in Conservation Priorities
Abraham and others 2010) also results in high levels of EFT rarity. This ecotone shows singular assemblages of species from tropical and nontropical biota (Gonzá lez-Abraham and others 2010) and a high diversity of distinctive lifeforms (Webb and Turner 2015). Finally, the surroundings of wetlands and mangroves in the Desert Region also showed rare EFTs, and both Mediterranean-type ecosystems and ecotones around wetlands are known to contain singular EFTs in other parts of the world (Cabello and others 2013). The lowest EFT rarity value was measured for the Tropical Region and southern desert mountains (Giganta Ranges), where heterogeneity and singularity are only introduced by the presence of endemism-rich evergreen pine forests at the highest altitudes (Leó n de la Luz and Domínguez-Cadena 1989).
As found at the species level (Riemann and Exequiel 2007;Lamoreux and others 2006), EFT richness and rarity were only correlated with a degree but did not always coincide in the peninsula. Such spatial aggregation between areas with both high EFT richness and rarity highlights their importance for heterogeneity and singularity in primary production.

EFTs for Setting Geographic Conservation Priorities
Three main conclusions can be drawn from our congruence analysis of the three approaches. First, our results highlight the importance of congruence areas as probable aggregated hotspots for all dimensions and scales of biodiversity, including diversity in essential ecosystem functions such as primary production dynamics. Areas with congruence reinforce their ecological and conservation value for the expansion of protected area networks (Lamoreux and others 2006). For instance, consistently prioritized areas of the Mediterranean mountains have been historically identified as a conservation gap based on plant diversity and endemism (for example, Garcillá n and Ezcurra 2003; Riemann and Ezcurra 2005). This congruence of the Mediterranean Region in North America suggests that some global biodiversity hotspots stand out not only as hotspots of endemism but also as heterogeneous and singular areas of ecosystem function, even if their identification does not consider ecosystem processes (Myers and others 2000). Second, our results indicate that traditional approaches may not identify all important areas of ecosystem functions (Meyer 1997) and may tend to better prioritize areas with rarity than those with richness in EFTs. Such an incidental focus of tra-ditional approaches on rare EFTs could derive from the dominant role that endemicity, often related to singular conditions, plays in conservation planning (for example, Myers and others 2000). It is interesting that heterogeneity in ecosystem functions has played a minor role (Lovett and others 2005) despite habitat heterogeneity fostering species adaptation and persistence (Hanson and others 2020). Third, our results also suggest that species diversity, as in hotspots of the Tropical Region mountains Ezcurra 2005, 2007), is not necessarily associated with rare or spatially heterogeneous ecosystem functions. In such areas, not high environmental heterogeneity but a long history of evolutive isolation under stable conditions has mainly driven speciation (Sundaram and others 2019).
Conservation efforts must employ spatially explicit and parsimonious ways to incorporate heterogeneity in ecosystem functions (Turner and Chapin 2005) to develop theories and tools that complement traditional planning and management actions (Possingham and others 2005). Our study shows how satellite-derived EFAs and EFTs of a focal ecosystem function (here primary production) offer tangible and biologically meaningful qualities of ecosystem functional heterogeneity (here EFT richness and rarity) that can complement traditional geographic priority approaches. EFAs and EFTs of focal ecosystem functions have already been used to assess the comprehensiveness and representativeness of protected areas others 2012, 2013) and of environmental observatory networks (for example, LTER, NEON, Ameriflux, and Mexflux; Villarreal and others 2018). Previous studies have also shown how EFAs and EFTs could facilitate conservation by capturing heterogeneity in the amount and timing of key ecosystem functions to model species distributions (for example, Tuanmu and Jetz 2015; Alcaraz-Segura and others 2017; Arenas-Castro and others 2018) and abundances (Arenas-Castro and others 2019) as well as provisioning, regulating, and cultural ecosystem services (Vaz and others 2020).

Caveats and Avenues for Future Research
The use of the EFT concept in geographic conservation is still subject to challenges. First, our satellite-derived EFT map characterizes the spatial heterogeneity of primary production dynamics. However, EFTs can also be identified from other remote sensing indices (for example, Ferná ndez and others 2010) to characterize the spatiotemporal heterogeneity of multiple ecosystem processes and functions at different scales to guide biodiversity and ecosystem services policies (Pettorelli and others 2018). Second, as the environmental observatory network expands, EFTs could be parameterized (for example, Mü ller and others 2014) and validated using ground measurements (for example, eddy covariance estimates of net ecosystem exchange; Villarreal and others 2018). Third, EFT richness and rarity maps illustrate diversity and spatiotemporal heterogeneity in the occurrence of ecosystem functions, but additional landscape indices could also elucidate the spatial arrangement (Fahrig and Nuttle 2005), connectivity, and lateral transfers (sensu Turner and Gardner 2015) of energy and matter fluxes at the landscape level. Fourth, our study does not assess interannual changes in EFAs, EFTs, or EFT richness and rarity, which could help reveal areas suffering from functional diversity homogenization, which is a planetary boundary that still needs evaluation (Steffen and others 2015). Fifth, the effects of spatial scale (grain and extent) on richness, rarity, and congruence with other biodiversity facets should be evaluated. Grain or cell size affects the magnitude, location, and spatial congruence of hotspots of species richness and endemicity (Rahbek 2005; Arponen and others 2012; McKerrow and others 2018; Daru and others 2020). The extent of the area under analysis may show that species-based priorities at one scale (for example, global) may or not overlap with those of other scales (for example, national or regional) (known as the parochialism effect; Pouzols and others 2014). EFT richness, rarity, and priorities depend on the extent considered but seem to be robust against sliding window sizes and the number of EFT classes defined (Appendixes 5 and 6). Future works should explore the effect of image pixel size (for example, with Sentinel-2 at 10 m/pixel), hierarchy in EFT classifications, and parochialism on the EFTbased approach. Finally, to test their effectiveness as ecosystem agnostic essential biodiversity variable candidates, EFT richness, rarity, and derived priorities should be compared to robust systematic conservation-based approaches that consider multiple facets of biodiversity, that is, compositional, structural, functional, and phylogenetic, in other ecoregions of the world (Pettorelli and others 2016).
In conclusion, the remotely sensed EFT approach can be used to incorporate the heterogeneity and singularity of ecosystem functions into geographic conservation priorities. Such an approach can support decision-making by offering supplementary arguments for the holistic conservation of biodiversity through the identification of key areas for multiple biodiversity facets (for example, the Mediterranean Region of Baja California) and of other areas important for ecosystem function that complement existing protected area networks (for example, mountainsides, and piedmonts with riverine systems in the Desert Region). Priority assessments based on essential variables related to ecosystem function cannot replace the use of very valuable systematic conservation approaches based on field records of species distributions to assess biodiversity status and change (Pereira and others 2013). However, our approach is useful to complement traditional priority setting, because is simple and based on only three satellite-derived meaningful descriptors of ecosystem functioning, facilitating computation and interpretation by managers and policymakers (Palumbo and others 2017). Future conceptual and empirical development and applications of EFTs should include other ecosystem functions, field validation, temporal changes in EFT diversity, and further metrics of heterogeneity across scales.