Testing drivers of acoustic divergence in cicadas (Cicadidae: Tettigettalna)

Divergence in acoustic signals may have a crucial role in the speciation process of animals that rely on sound for intra‐specific recognition and mate attraction. The acoustic adaptation hypothesis (AAH) postulates that signals should diverge according to the physical properties of the signalling environment. To be efficient, signals should maximize transmission and decrease degradation. To test which drivers of divergence exert the most influence in a speciose group of insects, we used a phylogenetic approach to the evolution of acoustic signals in the cicada genus Tettigettalna, investigating the relationship between acoustic traits (and their mode of evolution) and body size, climate and micro‐/macro‐habitat usage. Different traits showed different evolutionary paths. While acoustic divergence was generally independent of phylogenetic history, some temporal variables’ divergence was associated with genetic drift. We found support for ecological adaptation at the temporal but not the spectral level. Temporal patterns are correlated with micro‐ and macro‐habitat usage and temperature stochasticity in ways that run against the AAH predictions, degrading signals more easily. These traits are likely to have evolved as an anti‐predator strategy in conspicuous environments and low‐density populations. Our results support a role of ecological selection, not excluding a likely role of sexual selection in the evolution of Tettigettalna calling songs, which should be further investigated in an integrative approach.

Acoustic signals are subject to different evolutionary pressures.
Sexual selection, environmental conditions, morphological constraints and the presence of predators and/or parasitoids are just some of the factors that may act concomitantly or in opposition to each other (Endler, 1992;Wilkins et al., 2013). As a result of its influence on the divergence of sexual traits-which can lead to reproductive isolation-sexual selection has been regarded as an important engine of speciation (e.g. Butlin et al., 2012;West-Eberhard, 1983).
However, recent studies show that the combination of both natural and sexual selection is particularly powerful to initiate and complete speciation (Oneal & Knowles, 2013; van Doorn et al., 2009). For instance, song frequencies differ with body size, since in many animals there is a negative correlation between body length and dominant song frequency (Bradbury & Vehrencamp, 1998). This allometric relationship has also been demonstrated for cicadas, with sound frequencies seemingly being constrained by body size (Bennet-Clark & Young, 1994). Hence, ecological selection on body size may promote acoustic divergence and vice versa (Wilkins et al., 2013).
The acoustic adaptation hypothesis (AAH) postulates that acoustic signals diverge according to environmental physical properties that influence the attenuation and degradation of a signal during transmission through the atmosphere to maximize signal transmission and decrease signal degradation (Wiley & Richards, 1978. Consequently, animal sounds appear to adapt to a specific habitat in both frequency and temporal components (e.g. Badyaev & Leaf, 1997;Morton, 1975). Microclimate and vegetation structure are important in sound degradation due to atmospheric absorption and scattering, respectively. Sound attenuation increases with increasing temperature and decreasing humidity, while high-frequency sounds attenuate faster and are more scattered by foliage than low-frequency sounds (Morton, 1975;Wiley & Richards, 1982). Animals in hotter and drier microclimates should therefore produce lower frequency sounds than animals in more temperate or humid climates. Signals of wood warblers and bats were found to diverge along microclimate gradients due to climate influence in atmospheric absorption (Snell-Rood, 2012). Divergence was primarily at the spectral level, but bats also changed signal duration between seasons. Temporal and spectral attributes vary with vegetation structure because scattering increases positively with habitat closeness. Consequently, songs in closed habitats are predicted to have a longer duration, lower repetition rate and lower frequencies, with shorter notes and longer intervals, than songs in open habitats (Badyaev & Leaf, 1997;Ey & Fischer, 2009). While acoustic adaptation has been substantially studied in vertebrates (mammals, birds and anurans) with inconsistent results (review in Ey & Fischer, 2009), studies with insects are scarce and lend little support (e.g. Couldridge & van Staaden, 2004).
To investigate acoustic adaptation hypotheses in insects, one should control for phylogenetic history, and body size, and use populations of a species or closely related species occupying different habitats (Balakrishnan, 2016), which in most cases has not been available.
Males of Tettigettalna produce species-specific calling songs to attract females and are morphologically very similar. Sound is the main character used for species identification and is thought to play an important role in their reproductive isolation preventing hybridization among sibling species (Simões & Quartau, 2006).
Calling songs have a simple acoustic construction with no harmonic structure. Echemes are the main units of the calling song, grouping in species-specific rhythm arrangements that constitute phrases (Puissant & Sueur, 2010). Tettigettalna species exist only in the West Mediterranean Basin and occupy a variety of habitats. Most species are restricted to the Iberian Peninsula (see Figure 1) and are likely to represent a recent radiation (Costa, 2017). The current Iberian climate is very heterogeneous, influenced by both the Atlantic Ocean and the Mediterranean Sea, with North-South and East-West precipitation gradients, with increased precipitation towards North and West, while temperatures increase along a North-South gradient (Cunha et al., 2011;Rivas-Martínez, 2005). Following the Köppen-Geiger climatic classification, the peninsula is divided into two major climatic regions: (1) temperate with dry or temperate summer in the northeast, the west coast of Portugal, and the mountain regions and (2) temperate with dry or hot summers, covering most of the southern central plateau region and the Mediterranean coastal regions (Cunha et al., 2011). Additionally, local orographic and geological conditions allow for the formation of microclimates, from hot desert climate (e.g. Almeria, Murcia) to cold with temperate and dry summer climate found in the highest altitudes of the Sierra Nevada (Cunha et al., 2011). Tettigettalna argentata is the only species with a wide distribution (Gogala & Gogala, 1999;Hertach, 2008;Puissant & Sueur, 2010). This species is geographically structured into three or four genetically distinctive lineages (Costa, 2017;Nunes et al., 2014).
The remaining species occupy the south of the peninsula in a parapatric, patchwork-like pattern, except T. estrellae, whose distribution is limited to the North of Portugal (Sueur et al., 2004), and T. afroamissa which is endemic to the north of Morocco in the Rif and Middle Atlas Mountains . Nevertheless, several species' distribution ranges overlap.
Most species do not show host specialization and are observed in a broad diversity of plants within favourable habitats (Puissant & Sueur, 2010); however, some species are more often found in closed, arboreal vegetation while others almost exclusively inhabit open, shrubby or herbaceous vegetation types. The extent to which vegetation structure (open vs. closed) and local climate are of importance to the acoustic diversification of these species remains to be discovered.
Here, we explore the relationship between calling song variation and ecological divergence in Tettigettalna species, while accounting for body size and phylogeny, through the employment of phylogenetic comparative methods (PCM). These methods are an analytical toolkit used to study species in a historical framework to understand the mechanism behind their diversification (Garamszegi, 2014). Two types of data are required (1) an accurate phylogeny of the taxon under study and (2) interspecific phenotypic data for the same taxon (Garamszegi, 2014 High-phylogenetic signal has been construed as evolutionary conservatism (e.g. Ossi & Kamilar, 2006) or gradual evolution through time (e.g. Brownian motion) and is commonly connected with genetic drift, stabilizing selection or physiological constraints (e.g. Blomberg & Garland, 2002;Wiens et al., 2010). Low phylogenetic signal has been interpreted as evolutionary 'lability ' (e.g. Blomberg et al., 2003) or high evolutionary rates and is usually associated with adaptive radiations, divergent selection, or convergent evolution.
Inferring the evolutionary process from phylogenetic signal estimation, however, may be tricky and prone to misjudgements (Revell et al., 2008). Nevertheless, in conjunction with other analyses, it can provide insights into the underlying processes driving a particular trait evolution.
Our main goal is to investigate the relative role of stochastic (e.g.
genetic drift) and deterministic (e.g. ecological selection) processes in the divergence and evolution of the acoustic signals produced by 12 species of cicadas in a variety of microclimates and micro- Acoustic recordings were initially visualized and inspected for quality in Audacity version 2.1.0 (Audacity Team, 1999-2015. For each record, individual songs were identified according to the known description of each species' calling song (Costa, 2017;Puissant & Sueur, 2010;Sueur et al., 2004). Each individual song was then saved as a separate file. When necessary, the noise between echemes that could potentially interfere with echeme detection was manually removed. Recordings with background noise that could not be separated from the cicada sounds were discarded. Each song was then analysed with the R package warbleR version 1.

| Phylogeny
We obtained the most recent molecular phylogeny based on the concatenation method from Costa, 2017. This phylogeny includes all known species and identifies some genetically divergent lineages, which we used as species in our analyses ( Figure 2

| Phylogenetic signal
We tested the presence of phylogenetic signal in the morphometric, ecological and acoustic variables. For the continuous variables (all acoustic and morphometric variables), we estimated the evolutionary parameter Pagel's lambda (λ;Pagel, 1999), using the phylosig function in the R package phytools version 0.7.47 (Revell, 2012). We chose Pagels's lambda over other indices as it proved to be the least affected by the number of species in the phylogeny (Münkemüller et al., 2012). Pagel's λ follows a Brownian motion (BM) model of trait evolution, which assumes that phenotype traits evolve through a series of continuous, random, independent steps and variance between species accumulates in proportion to shared phylogenetic history (Felsenstein, 1985). λ normally varies between 0 (no phylogenetic signal) and 1 (high phylogenetic signal).
A value of lambda equal to 0 indicates that traits covary less than expected under BM (i.e. independently of phylogenetic history), whereas a value of lambda equal to 1 indicates that traits covary as expected under BM (i.e. closely related species are more similar to each other).
For the binary habitat variables (calling site and NDVI), we esti-

| Correlated evolution of calling song
The evolutionary association between each acoustic variable and the morphological and ecological traits was analysed using phylogenetic generalized least squares analyses (PGLS; Martins & Hansen, 1997 Clade 2 phylogenetic covariance matrix in the error structure of the model. The expected covariance matrix is proportional to the amount of shared evolutionary history between the species in a given phylogeny (assuming a BM model of evolution) and likelihood methods can be applied to find transformations of the matrix that improve the fit of the model. The most used transformation is the λ transformation. When λ = 1, the model assumes a BM model of evolution, while values of λ < 1 indicate that the covariation between the traits is lower than expected under Brownian evolution (Symonds & Blomberg, 2014). We used the function pgls in the R package caper with λ transformation to test each hypothesis of correlated evolution, separately. Temporal acoustic variables and pronotum length were log-transformed, and dominant frequency was squaretransformed to meet the statistical requirements of normality and homogeneity of the model's residuals (Freckleton, 2009). Bioclimatic variables were not transformed. For each model, we inspected the distribution of the residuals to assess their normality and homogeneity (Freckleton, 2009).

| Mode of acoustic evolution
To visualize evolutionary patterns in multivariate phenotypic (acoustic) space, we plotted a 'phylomorphospace', which superimposes a phylogenetic tree into a subspace defined by the principal components of phenotypic variation. Phylomorphospace is a common tool used in the study of macroevolutionary processes in morphological change (Adams & Collyer, 2019) and here we extended its use to acoustic space. We used the phylomorphospace function in the R package phytools to project the phylogeny into acoustic space, ancestral states of the internal nodes were estimated by maximum likelihood. Since dominant frequency did not show to be important between species and for an easier interpretation, only the temporal acoustic variables were used (see section 3).
To discern the evolutionary processes that might influence the diversification of Tettigettalna acoustic signals, we assessed the fit of seven evolutionary models to our acoustic variables. We fitted three single-rate models and four multi-regime models to our data. For the single-rate models, we considered: Brownian-Motion (BM), singleoptimum Ornstein-Uhlenbeck (OU) and Early-Burst (EB) models. The BM model (Felsenstein, 1985) describes trait evolution as a random walk wherein the covariance structure is defined by the phylogeny (i.e. closely related species are more similar) and can be described with two parameters, the evolutionary rate parameter (σ 2 ) and the mean starting value of the trait (z(0)). This model is commonly associated with evolution by neutral drift; however, a BM pattern of evolution can arise by other processes such as selection in a quickly fluctuating environment (Hansen & Martins, 1996). The OU model (Butler & King, 2004;Hansen, 1997)  variance, all models were fitted with the standard error calculated from each species' data. All single-rate models were fitted using the function fitContinuous in the R package geiger (Pennel et al., 2014) and the multi-regime models with the function ouwie from the package OUwie version 2.2. (Beaulieu & O'Meara, 2020). Model fit was assessed using Akaike's Information Criterion corrected for small sample sizes (AICc) and Akaike weights (AICw). An ΔAICc threshold of 4 was used to infer strong support (Burnham & Anderson, 2002).

| Sampling and habitat characterization
Our final data set contained 1558 songs from 255 individuals spread across 14 species (Table 1). As expected, most species were found calling from a variety of substrates. Nine species called predominantly from trees and five from shrubs ( Figure 3).

| Phylogenetic signal
Estimates of λ showed that only one acoustic variable, echeme rate, had a strong phylogenetic signal significantly different from 0 (p < 0.05;

| Body size
Phylogenetic generalized least squares analyses analysis revealed a significant negative correlation between body size (PL) and dominant frequency (PGLS: R 2 = 0.25; p = 0.04; Table 3; Figure S1A). There seems to be a great influence of T. josei in this relationship ( Figure S1A).
There was no association between any temporal acoustic variables and body size ( dominant frequencies, and that body size does not influence temporal patterns, which agrees with our expectations that frequencies are constrained by the size of the sound production systems.

| Climate
Nine bioclimatic variables were retained for analysis (Table 4). We found a significant relationship between two temporal acoustic variables (number of echemes and interval duration) and some of these (

| Habitat structure
Phylogenetic generalized least squares analyses analysis between acoustic variables and habitat openness/closeness (NDVI , Table 6) revealed that the number of echemes (R 2 = 0.56; p = 0.001) and echeme duration (R 2 = 0.26; p = 0.036) are positively correlated with NDVI. Figure 4a,b show that species in low NDVI habitats have a lower number of echemes and longer echeme duration. The exception was T. estrellae and T.h. galantei Type II that showed a pattern similar to low NDVI species. We also found a correlation between calling site (Table 7) and interval duration (R 2 = 0.32; p = 0.021).
Shrubby species have a higher interval duration than arboreal species ( Figure 4c). The exception to this pattern is T. josei, with an opposite trend. Overall, species in closed habitats (high NDVI or ar-

| DISCUSS ION
Our approach in this paper was to investigate the relationship between calling song variation and ecological divergence in a

TA B L E 6
Phylogenetic generalized least squares model of the relationship between each acoustic variable and NDVI (high/low). phylogenetic framework while accounting for body size. We attempted to infer the evolutionary history of six acoustic parameters, and if the major forces shaping their divergence were stochastic (e.g.

F I G U R E 4 Correlated evolution between calling site and NDVI
genetic drift) or deterministic (e.g. ecological adaptation) in nature.
We included in our analysis predictions from the AAH, which postulates that acoustic communication should adapt to the physical properties of the environment in which it is produced in a way that maximizes transmission and minimizes degradation. If stochastic processes are the main driver of acoustic divergence, we expect to find no correlation between acoustic and ecological divergence and a significant phylogenetical signal; if on the other hand ecological adaptation was the most influential mechanism then acoustic divergence should relate to ecological variation and not be phylogenetically structured. Our results indicated an interplay between stochastic and deterministic processes and different evolutionary histories for each acoustic variable. We found little support for the Acoustic Adaptation Hypothesis.

| Phylogenetic signal and mode of evolution
We found a statistically significant phylogenetic signal in one temporal acoustic variable (echeme rate) and one ecological variable (calling site). The remaining acoustic variables, body size and the ecological variable NDVI did not display a significant phylogenetic signal. These results are in partial disagreement with the literature. Fonseca et al. (2008) contrasted the phylogenetic history of nine cicada species (including four Tettigettalnas) with their calling song patterns and the morphology of the song production apparatus.
The authors found that the spectral traits of the song, which are strongly dependent on body size (Bennet-Clark & Young, 1994) and the mechanisms of the tymbal (Fonseca & Popov, 1994) showed a strong correlation with phylogenetic relationships. On the contrary, the calling song temporal pattern is controlled by a highly plastic nervous system, which consequently exerts fewer biomechanical constraints on song evolution. As a result, the temporal pattern was independent of the phylogenetic history (Fonseca et al., 2008). A reasonable explanation for the lack of signal in the spectral variable is that our data are comprised of a single genus. Had we included other genera such as Cicada or Tibicina, both dominant frequency, and body size would probably present some degree of phylogenetic signal, as shown by Fonseca et al. (2008). Indeed, we found the expected correlation between body size and dominant frequency; around 25% of the variation in dominant frequency is explained by body size when taking phylogeny into account. Frequency ranges seem to be genus-specific and constrained by body size (Fonseca et al., 2008), being very similar between Tettigettalna species, except for T. josei, the smallest and earliest to diverge in the genus, which calls at much higher frequencies than predicted for its size.
Moreover, results from the model fitting analysis suggest an earlyburst model of evolution for the dominant frequency variable. It is unlikely that dominant frequency is used by female and male cicadas to distinguish conspecifics from heterospecifics of similar body size.
Therefore, they probably rely on temporal parameters when found in sympatry. How females perceive these acoustic cues and whether they have an impact on mate preferences remains to be tested.
The discovery of a significant phylogenetic signal on echeme rate indicates that closely related species have similar values for this variable. A high phylogenetic signal is usually synonymized with evolution by neutral genetic drift. However, given that similar phylogenetic signals can be produced by several different evolutionary processes (Revell et all., 2008), results should be interpreted with caution and in conjunction with other analyses. The PGLS analysis did not reveal any correlation between echeme rate and the ecological or morphological variables. Furthermore, the model fitting analysis did not deliver strong support for a single model, but the BM and the two BMS models were similarly informative. The overall results point to evolution by genetic drift in this variable. The variable call duration had no phylogenetic signal; no correlation with any ecological or morphological variable and the model fitting analysis indicated a Brownian Motion model of evolution with strong support. Given these results, the absence of a phylogenetic signal should not in this case be interpreted as indicative of a deterministic process. Revell et all (2008), explicitly warn that a low phylogenetic signal can be produced for all but one (constant-rate genetic drift) of the processes tested in their study. They found that a decreased phylogenetic signal can be produced under genetic drift when the rate of evolution was initially low but increased over time since this process tends to increase variation across the tips of the tree without a corresponding increase in the covariance among taxa. Consequently, the evolution of this variable seems once again best explained by genetic drift.
Three variables (number of echemes, echeme duration and interval duration) displayed a correlation with some ecological variables and no phylogenetic signal, suggesting ecological selection.  (Landis & Schraiber, 2017). Fluctuating natural selection can also produce a low phylogenetic signal when the rate of fluctuation is low (Revel et al., 2008).  (Grace & Shaw, 2004) and field crickets (Walker, 2000) with higher developmental temperatures resulting in faster songs.

| Ecological selection and the AAH
However, the sound production systems of the Orthoptera (stridulation) and the Cicadidae (by means of a tymbal mechanism) are very different (Claridge, 1985;Forrest, 1982), as well as the larval stage duration, that in cicadas can be up to 17 years underground (Williams & Simon, 1995). The most likely explanation is that the correlation found between climatic variables (BIO 2,4 & 5) and some acoustic variables (number of echemes & interval duration) is a by-product of the correlation between climate and vegetation.

It is well documented that climatic conditions influence vegetation
and plant distribution (e.g. Gavilán, 2005). Temperature seasonality (BIO4), the bioclimatic variable with greater influence on the acoustic structure, is a measure of continentality. In the Iberian Peninsula, it decreases circularly from the interior to the coasts, being particularly low on the Atlantic shores and high in the centre and south of Spain (Andaluzia) and its mountainous regions. The latter regions, particularly the high altitudes of the Sierra Nevada and Sierra de la Contraviesa are also regions of low NDVI, dominated by shrubland. A recent study using climatic data to seek the drivers of floristic composition differentiation in the thermophilous deciduous oak forests in Eastern Europe found that temperature seasonality was a major variable separating several clusters of forest types (Goncharenko et al., 2020). Indeed, we found an association be-  (Sueur & Aubin, 2002). Playback experiments with C. orni showed that interval duration may be essential for species discrimination (Simões & Quartau, 2006). The acoustic adaptation hypothesis refers especially to long-distance signals (Ryan & Kime, 2003), since these are more susceptible to degradation. Tettigettalna frequencies are genus-specific, while temporal parameters are species-specific (Fonseca et al., 2008). It seems likely that these species use frequency for long-distance communication and temporal proprieties for short-distance communication, thus circumventing signal degradation. Short-distance signals are not under strong natural selection and can even be selected for fast degradation to avoid predators and parasitoids (Endler, 1992).
In fact, is well documented that predation is an important pressure in the evolution of acoustic and visual signals. Predation of cicadas has been reported in birds (Pons, 2020;Williams et al., 1993) small mammals (Krohne et al., 1991), lizards (Han et al., 2021) and several arthropods such as spiders (Suzuki & Mukaimine, 2021) and ants (Whitford & Jackson, 2007). Cicadas are also potential hosts for phonotactic parasitoid flies. Emblemasoma species (Diptera: Sarcophagidae) were found parasitizing different cicada species such as Tibicena pruinosa and T. chloromera (Farris et al., 2008), T.
dorsatus (Stucky, 2015) and Okanagana rimosa (Schniederkötter & Lakes-Harlan, 2004 Population density is an important factor in predation avoidance. Williams et al. (1993) found that avian predation in Magicicada is inversely density-dependent due to predation satiation, reaching 100% for individuals emerging too early or too late in the season.
Periodical cicadas are much easier to catch (both by humans and birds) than non-periodic cicadas, having been labelled as 'predatorfoolhardy' (Lloyd & Dybas, 1966). The synchronized life cycle and large aggregations produced by these periodical cicadas are thought to have evolved precisely as an anti-predation mechanism (Karban, 1982;Lloyd & Dybas, 1966 Signal evolution is a delicate equilibrium between local predation risk and sexual selection (Endler, 1992). Where the predation risk is higher and the strength of sexual selection lower, signals evolve to be more cryptic (Endler, 1992) or even degrade faster at shorter distances (Wiley & Richards, 1982). An extreme example is the case of the Pacific field cricket, Teleogryllus oceanicus, where a 'flatwing' mutation that causes sexual signal loss, spread in <20 generations to 5090% of the population in two Hawaiian Islands (Zuk et al., 2006(Zuk et al., , 2018. The 'flatwing' mutation is a rapid adaptive response to a deadly acoustically orienting parasitoid fly, whose range overlaps with T. oceanicus only in these islands. 'Flatwing' males showed increased phonotaxis to normal-wing males and acted as 'satellites' to circumvent the difficulty of attracting females without song. A pre-existing plastic behaviour coupled with strong natural selection against calling males allowed this otherwise maladaptive mutation to become established (Zuk et al., 2006) Sexual selection has been regarded as a major driver of speciation (e.g. Fisher, 1930;Lande, 1981Lande, , 1982, this is especially so, for sexual traits, as is the case with acoustic signals (e.g. review Wilkins et al., 2013). More recently, it has been argued that speciation by (pure) sexual selection is highly unlikely and that sexual selection contributes more effectively to speciation alongside other processes, like ecological selection (Ritchie, 2007). Cornwallis and Uller (2009) argue for an 'evolutionary ecology of sexual traits' with a focus on how the interplay between environmental heterogeneity and phenotypic plasticity influences the evolution of sexual phenotypes.

| CON CLUS ION
Acoustic signals play an important role in species divergence and speciation. As such the study of their evolution is of great importance to understand the evolution of the species that rely on them for mate attraction and communication. Discerning the evolutionary process behind acoustic divergence is not always easy or even possible. By their very nature, acoustic signals are subject to a myriad of different pressures, which may act in isolation, in synergy or antagonistically to each other. Here, we present the most extensive test for acoustic ecological adaption in insects using a phylogenetic framework, including predictions of the AAH. Our results show some support for ecological adaptation at the temporal but not the spectral level, with no support for AAH.
Acoustic divergence in this group of cicadas seems to have been influenced by both stochastic and deterministic processes, acting distinctively in different variables. The most influential explanatory variables-BIO4, NDVI, and calling site-seem to reflect some measure of habitat openness. Knowing which one is the proximal variable that directly affects song structure, and which are the distal variables that act through correlation with the proximal, requires further testing. However, all three influence song structure in the same way: Songs in open habitats are defined by a low number of echemes and long echemes and intervals. These characteristics are contrary to what is expected under the AAH and can cause them to degrade more easily in these habitats.
These traits are likely to have evolved as an anti-predator strategy in conspicuous environments and low-density populations.
We hypothesize that the interplay of genetic drift (for most of the phylogenetic history) and local ecological adaptation (when faced with high-risk habitats) played an important role in acoustic divergence in this genus.
Considering these results, it becomes relevant to verify whehter these patterns can be generalized to a higher taxonomic level. For

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.14133.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at: https://doi.org/10.5061/dryad.jh9w0 vtfk.