Quantifying phenological diversity: a framework based on Hill numbers theory

Phenology as a continuous variable: wavelet time series analysis

Wavelet transform detects the frequency spectrum of discrete time series data and fits a smoothed curve to them (Schmidt & Skidmore, 2004; Sifuzzaman, Islam & Ali, 2009). This analysis is based on the comparison of the similarity between a scaling function (which can be stretched, shrunk, and shifted in time) and the original time series (Torrence & Compo, 1998). With these comparisons, a matrix is constructed that contains the fits of the scaling function to the time series, in which the total sum of each column in the matrix produces the smoothed curve. In the case of phenological data, we can use such a smooth curve as an approximation to a continuous phenological pattern, in which the abrupt changes caused by sampling effort and protocol are smoothed out, in line with the assumption that phenophases start and finish gradually rather than abruptly (Fig. 1).

Wavelet transform analyses require specifying two parameters. The first parameter is the scaling function. The most common scaling functions are Morlet, Paul, DOG, Biorthogonal and Mexican hat, and these differ mainly from each other in the signal resolution calculation (González-Nuevo et al., 2006; Singh, Singh & Sharma, 2011). For our purposes, the Morlet scaling function is preferable, as it is optimal for data that cannot be directly interpreted, and time series with unknown frequencies and scales (Percival & Walden, 2000). The second parameter is the attenuation threshold (τ), which determines the phenological curve steepness. Values close to 0 translate into a highly smoothed curve; conversely, values approaching infinity translate into a wiggly curve, similar to the raw time series data. A value of 2 is used as a standard for wavelet analysis in mathematics software (e.g., Matlab and wavScalogram R package). Wavelet analysis must be performed on the discrete phenological curve or time series of each studied species in the community; thus, the number of total smoothed phenological curves that must be calculated equals the number of species in the community (S).

Quantifying Phenological Hill numbers (PD)

Our phenological diversity measure is based on the Hill numbers-based attribute diversity framework (Jost, 2007; Chiu & Chao, 2014; Chao, Chiu & Jost, 2014). Hill numbers are a set of metrics that have two major advantages over other diversity indices: (1) the interpretation of the diversity values are consistently the same, and (2) the sensitivity regarding abundant and rare species or traits can be regulated with a parameter (q). This q parameter directly determines the sensitivity of the diversity measure (^qD) to the relative abundances of species occurring in the community. Although q can take any non-negative real number, ecologists typically consider three values: 0, 1 and 2 (Chao, Chiu & Jost, 2014). When q = 0, the relative abundance of species is overlooked, and the measure simply represents S (i.e., ⁰D = S). When q = 1, S is weighted by the proportions of the species abundances and can be interpreted as the effective number of species equally abundant within an assemblage, which is equivalent to exp(H) (i.e., the exponent of Shannon’s entropy; Jost, 2006). Finally, when q = 2, the diversity values favor the most abundant species, and the less abundant or rare species are almost not accounted for; consequently, ²D can be roughly interpreted as the effective number of dominant or the most abundant species in the community (Jost, 2006). In addition, Hill numbers are consistent with basic diversity concepts like evenness and dominance (Chiu & Chao, 2014; Chao, Chiu & Jost, 2014). Additionally, Hill numbers are expressed in intuitive units of effective number of species, and they can be directly compared across orders of q to gain information on the dominance, community traits and comparisons among different species assemblages. Finally, Hill numbers theory can be generalized to taxonomic, phylogenetic, and functional diversities (Chiu & Chao, 2014; Chao, Chiu & Jost, 2014; Chao et al., 2019). Here, we use this framework to measure PD in a community.

Phenological diversity assessment through the species-pairwise distance framework

The measure that we present here is based on a pairwise overlapping distance, following the same logic by Chao et al. (2014) in developing their functional diversity measure under the assumption that each species has specific phenological curves (Chao, Chiu & Jost, 2014; Chao et al., 2019). The distance we used is based on the Morisita-Horn index modified to measure the amount of overlap between pairs of phenological curves (Luna-Nieves et al., 2022). Let O_ij be the pairwise overlapping distance between the continuous phenological curves of the i-th and j-th species, defined as (1)${\displaystyle o_{ij}=1-\frac{2\int z_i\left(t\right)z_j\left(t\right)dt}{\int z_i{\left(t\right)}^{2}dt+\int z_j{\left(t\right)}^{2}dt}}$ where z_i and z_j are the smoothed continuous-over-time phenological curves of species i and j, respectively and the integral is calculated over the studied time interval. O_ij ranges in the [0, 1] interval, with O_ij = 1 when curves fully overlap (Fig. 2C), and O_ij = 0 when curves show no overlap (Fig. 2A). When completely overlapping in time (Fig. 2C), the species belong to the same phenological group; when they partially overlap, they partially belong to the same phenological group (Fig. 2B). A third case corresponds to the scenario in which the phenological curves do not overlap (Fig. 2A), which represents the existence of completely different phenological groups. Therefore, our approach to measuring phenological diversity is based on the pairwise overlap distance measured through the Morisita-Horn index (Horn, 1966; Garratt & Steinhorst, 1976).

Considering that both deterministic and stochastic variables contribute to the phenomena regulating phenology, similar phenological curves between pairs of species might suggest that they have similar environmental requirements, a common response to the same environmental cues, or similar evolutionary constraints. By contrast, dissimilar phenological curves suggest that the species have different environmental requirements or might simply reflect historical competitive displacement between species (Slagsvold, 1977; Baselga, Gómez-Rodríguez & Lobo, 2012). As the phenological curves may vary among time in position and shape. The overlapping area of the phenophase curves between pairs of species over the entire time series can provide an indirect and prospective measure of the interactions between phenophases only in the time frame measured. Because our framework of phenological diversity is based on the approach of “attribute diversity” (Chao, Chiu & Jost, 2014; Chao et al., 2019), which is a robust extension of Hill numbers, it can be applied to measure species traits and their diversity in orders of q. Within this framework, phenological Hill numbers (PD) can thus be interpreted as the pairwise phenological distance between species (in units of equally abundant and distinct species with distinct phenological group) occurring in an assemblage within the time interval over which it was measured.

To construct our phenological diversity measure we need to consider a standardizing factor: sum of the relative overlap of phenological curves, denoted as Q. This factor is calculated as follows (Chao, Chiu & Jost, 2014; Chao et al., 2019): (2)${\displaystyle Q={\sum }_{i=1}^S{\sum }_{j=1}^S{O}_{ij}{p}_i{p}_j,}$

where O_ij is calculated as in Eq. (1), and p_i is the relative intensity or abundance of the phenological event measured on species i, defined as (3)${\displaystyle p_i=\frac{\int z_i\left(t\right)dt}{\sum _{j=1}^S\int z_j\left(t\right)dt}.}$ Finally, the phenological Hill numbers of order q, ^qPD, is calculated as (4)${\displaystyle \text{ˆqPD}={\left[{\sum }_{i=1}^S{\sum }_{j=1}^S\frac{o_{ij}}{Q}{\left(p_i{p}_j\right)}^q\right]}^{\frac{1}{2\ast \left(1-q\right)}}.}$ Given that when q = 1 the exponent $\frac{1}{2\ast \left(1-q\right)}$ is undefined, we decided to use the same approach as Chiu & Chao (2014). Therefore, (5)${\displaystyle \text{ˆ1PD}=exp\left(-\frac{1}{2}{\sum }_{i=1}^S{\sum }_{j=1}^S\frac{o_{ij}}{Q}{p}_i{p}_{j}log\left(p_i{p}_j\right)\right).}$

Simulations and field data

To illustrate the utility of ^qPD, we generated simulated data sets reflecting different community scenarios and applied our measure to these and one additional real data sets. Although these simulations are unlikely to occur in nature, they provide a glimpse into a real community phenological pattern (Fig. 3). Simulations vary in several community traits, such as: the degree of overlap, differential or even in intensity, and total number of species (i.e., S) (Fig. 3).

As previously explained, ^qPD values represent the “phenological Hill numbers” and thus quantify the diversity of different phenological curves in a given assemblage. The contribution of the phenological curve of each species is considered unique and equally distinct from each other; thus, ^qPD values always range from >0 (Fig. 3) to the total number of phenological curves measured (S, when there is no overlap at all between them). Therefore, phenological Hill numbers values have lower values when q increases.

In addition to the simulated data sets, we performed an analysis of ^qPD based on field data for the breeding phenological activity of an amphibian community from Madagascar (Fig. 4; S = 40; time period = 360 days; frequency = daily; Heinermann et al., 2015).

All analyses were performed in R v. 4.1.2 (R Core Team, 2021) using the DescTools (Signorell et al., 2020) and wavScalogram (Benítez, Bolós & Ramírez, 2010) packages. We provide the script for calculating the phenological diversity measure in the Supplementary Material.