Cicadas (Cicadoidea) are an interesting focal group for ecological studies. They are abundant in diverse natural and artificial habitats around the world, including in the Americas, Africa, Eurasia, Australia, New Guinea, Borneo, Madagascar, and most non-arctic large islands (Marshall et al., 2018). Juvenile cicadas or nymphs are root feeders, feeding on the xylem fluid of plant species for several years, and thus are potential indicators of the integrity of the root system, and of the water and nutrients that host plants are absorbing (Williams and Simon, 1995). Nymphs also need uncompacted substrates to construct their burrows and feeding cells; their abundance being indicative of well-structured and healthy soils (Boulard, 1965; Smith and Hasiotis, 2008). Consequently, cicada nymphs have been used so far as bioindicators of human impacts on the soil, including its compaction (Moriyama and Numata, 2015), degradation (Hou et al., 2018) and contamination by mercury (Zheng et al., 2010).

Estimated biomasses of North American periodical cicadas (i.e. those whose adults appear synchronically every 13 or 17 years) are amongst the highest of any animal in terrestrial ecosystems (Karban, 2014). Yet, non-periodical cicadas (i.e. those whose adults can appear each year) can also reach tens of thousands of individuals per hectare (Patterson et al., 1997). As a result, the growth of trees holding large cicada numbers is limited through the costs imposed by chronic root herbivory (Karban, 1980; Yang and Karban, 2009). However, plant growth can also benefit from fertilization due to carcass deposition and belowground water transport by high-density cicada populations (Andersen, 1994; Yang, 2004). The aggregate nature of these populations could lead to increased spatial variability in vegetation growth and recovery after disturbances.

Cicada adults are important prey for arthropods (Suzuki and Mukaimine, 2021) and vertebrates (Pons, 2020). Their abundance may enhance the survival and reproduction of predator species (Karban, 1982; Koenig and Liebhold, 2005) and thereby contribute to the regulation of insectivorous populations. Nonetheless, high cicada population densities, as found in emergence years of periodical cicadas, are able to satiate the local community of avian predators, thus increasing the likelihood of adult cicadas escaping predation (Karban, 1982). Furthermore, dead cicadas can provide a significant input of insect biomass into forest soils (Yang and Gratton, 2014). The specialised relationships of cicadas with fungi and arthropod parasites and bacterial endosymbionts have also been described recently (Liu et al., 2018; Simon et al., 2022).

In temperate latitudes, cicadas have a relatively short and predictable period of emergence from their belowground life (Smith et al., 2006), whereas in tropical ecosystems adults occur for longer periods of time (Lee et al., 2010), or even year round (Hart et al., 2015). When emerged nymphs moult into adults, they leave exuviae (skins) attached to the lowest vegetation strata. Exuviae mark the precise location where nymphs have occurred belowground for several years (Boulard, 1965), and can therefore be used to gather information about cicada juvenile life and environment and be accurate soil bioindicators. Ecological relevance, indicator use, accessibility of exuviae, population densities and global distribution, all make Cicadoidea a candidate focal group for promising ecological studies. However, a standardised method for quantifying cicada exuviae in natural habitats is lacking.

Insect moult is the casting of the external cuticle that functions as the body exoskeleton. The skin or shed left after a moult is known as the exuvia. It is composed of chitin (a polysaccharide that amounts to 20-50% of the cuticle dry weight), proteins, lipids and minerals (Chapman, 2012). Lasting for several weeks before degrading, insect exuviae can be used in ecological studies. Insect exuviae have been counted and collected by different methods, especially in dragonflies and damselflies (Odonata; Raebel et al., 2010, but see Bried et al., 2012) and midges (Diptera: Chironomidae; Ruse, 2010). Their sampling is not weather dependent, as it is for adults, it provides an accurate estimation of sex ratios and numbers of emerged imago (Dolný et al., 2014), and does not affect population numbers or the welfare of sampled animals. Methods used to collect dragonfly exuviae include unspecified sampling on different substrates or areas (Raebel et al., 2010; Dolný et al., 2018), line transects for which area but not time is specified (Dolný et al., 2014), and time-limited sampling of sites of an unspecified area (Hardersen et al., 2017). Most studies do not describe methods used in detail and standardization was lacking until the work of Arguel et al. (2022) for riverine dragonfly communities.

In cicadas, the last moult separates the final instar nymph, living belowground, from the adult living aboveground (Boulard, 1965). Final instar nymphal exuviae are found attached to grasses, shrubs and tree trunks, in well-exposed supports, since the newly emerged adult needs open space to unfold its wings and start the first fly (Boulard and Mondon, 1995). Body length of cicada exuviae (between the head and the last abdominal segment) oscillates from less than one to a few centimetres (Motta, 2003; Logan and Connolly, 2005; Lee et al., 2012). Direct sunlight and rainfall seem to accelerate their degradation, although this is poorly understood. Detectability of exuviae may vary according to their size, colour contrast with the background, density of low vegetation and whether they remain attached to the moulting support or have fallen from it (pers. obs.). Taxonomical and sexual identification of cicada exuviae by morphology or genetic analyses allows population studies to be conducted (Lee et al., 2012; Pons, 2015).

In contrast to the numerous studies based on odonata exuviae, cicada ecological studies have seldom used this source of information. More often, traps are used to capture newly emerged adults for estimating their feeding ecology (Callaham et al., 2000) or their density and emergence phenology (Smith et al., 2006; Whiles et al., 2019). Counting emergence holes or recording sound levels are also means to estimate cicada densities (Patterson et al., 1997; Maier, 2015). Cicada exuviae, specifically, have been counted in small quadrats of 0.125 m² along transects (Andersen, 1994), in 100 m² plots distributed on two-dimension grids (Lee et al., 2010) or in individual plots of variable size (35 to 14,000 m²) (Kim et al., 2014; Moriyama and Numata, 2015; Sato and Sato, 2015). Although several studies collect cicada exuviae, methodological limitations, such as poor sampling design, unknown sampled area or time, unmeasured detectability of exuviae, etc., usually prevent estimating population densities. To our knowledge, no cicada exuviae counting method exists that is described in detail and standardized so it can be used for monitoring cicada emergence in most regions and habitat types and estimate population densities.

We describe here a Cicada Monitoring protocol using Exuviae Transects, Cicada-MET. It allows abundances of emerging cicadas to be estimated in natural or semi-natural habitats. The method can be used to obtain long-term temporal trends, to monitor populations after natural or human induced disturbances, and to compare among experimental treatments, among habitats or, given some conditions, among taxa. This paper describes the different steps and gives recommendation for the use of the method and ideas for further studies. The quantitative aspects of the method rely on a 10-year database originated from the sampling of transects and plots in Catalonia (NE Spain).

The method was developed in the frame of a fire ecology study and was inspired in transects conducted previously (Pons, 2015). Data used to describe procedures and test assumptions were gathered from burned and unburned Mediterranean-type forests. Cicada-MET has been used in the field by trained observers who have conducted transect counts for ten years (PP and JMB) and plot counts for one year (AP and Gemma Banchs).

To conduct Cicada-MET, we recommend following seven consecutive steps that are summarized in Table 1 and described in detail in following sections. Additionally, a sample of transects can be counted twice successively to account for species detectability (see section 4.2). The procedure works if phenology is predictable and its interannual changes are moderate. It can be used for a focal species or for several species (with enough sample size in the study area), separately. It combines a time-consuming method (plots) that has to be applied once per study region (and ideally repeated each 5-10 years to account for longer term variability in phenology), with a rapid method (transects) that allows a high number of replicates each year.

For transects to be comparable, path characteristics have to be fixed in advance. This includes defining the range of track type, width and embankment occurrence that can be sampled. We avoided paved roads and paved margins since pavements can affect root systems. Embankments are often unsuitable for sampling because of high slopes and absence of cover. It is therefore better to avoid those track sections or to sample from behind the road banks.

Another aspect to consider is the human impact on roadside vegetation that may include clearing of grass and shrub layers, tree thinning, livestock grazing or controlled burning. When these disturbances take place before the start of the emergence season, they have limited effects on counts. However, if any of them takes place during the emergence season, the affected transect should not be counted (or the counts not pooled to the rest) because part of the existing exuviae would be lost.

The usual spatial sampling design choices apply to this method. These include an adequate spatial replication with interspersion of treatments (Hurlbert, 1984), and selecting a limited and well-known range of environmental characteristics for the study sites.

Regarding the temporal design for transect censuses, one visit per year may be sufficient if repeated monitoring of plots has been performed across an emergence season. The selection of sampling dates for plots can rely on expert knowledge of the emergence phenology of the species in the region, on scientific publications, on open biodiversity databases, or on previous fieldwork. Sampling dates for transects should start after the peak of emergence of the focal species and take place for no more than one month, although the extent of the sampling period should depend on local meteorological conditions. The percentage of exuviae fallen from their moult support can therefore be used as an indicator of the effect of rain and wind (but also other potential perturbations) on exuviae persistence (See under 4.3. Loss of exuviae due to weather).

The dates of counts have to be distributed among the different treatments/habitats so that temporal variation related to meteorology affects them equally. In practice, the optimal solution is to cover approximately the same number of transects per treatment or habitat to be compared, each sampling day.

Transects follow paths or tracks along natural or semi-natural vegetation (e.g. forests and tree plantations). Both sides of the path can be sampled by two observers simultaneously, or by the same person on the way out and on the way back. A standard area for transects is 200 m², with 100 m of length and 1 m of width on each side of the way. Transect width can be checked using a 1-m long stick while walking. We do not recommend sampling more than 1.5 m of width, because the visibility of exuviae rapidly decreases with increasing distance from the observer.

Observers need to be trained in finding exuviae ( Appendix A shows the exuviae morphology and length for the five species found in this study) in the vegetation and on the ground, before transects can be done reliably. Observers walk slowly surveying all visible vegetation and ground, with the help of the stick. The walking speed depends on the structure of the low vegetation and the percentage of bare ground, a usual speed is 0.4-0.6 km/h. However, the overall speed to cover a transect is lower, because it diminishes with the number of exuviae found and the number of field variables collected (see under 4.1. Sampling speed in transects). A 200-m² transect typically needs 30 to 90 minutes to be completed when exuviae are collected and three field variables (GPS position, plant species as support and height of exuviae) are measured. The time needed to complete the count is halved if two observers walk on each side of the path simultaneously.

Identifying cicada species through exuviae depends on coloration, size, and shape of the exuviae structures (Motta, 2003; Lee et al., 2012, Appendix A ). As a result, conducting a preliminary study becomes imperative in the absence of species identification keys. While some exuviae can be identified solely at the genus level, such as Cicadetta spp., specific determination is still feasible provided the cicada fauna in the study areas is well-documented. For example, if only one species from the genus is present, or if adults have been collected or their sound recorded near the transect. Moreover, a single exuvia provides a sufficient quantity and quality of DNA for genetic analyses and species identification (Bouwer et al., 2014). Additionally, distinguishing between males and females is feasible by examining the shape of the genitalia in the terminal abdominal segments of the exuvia (Motta, 2003).

It is important to clearly define transect limits. As it is difficult for markings to last for a long time on roadsides, we used a standard GPS device for an approximate location (usual precision +/- 3 to 7 m) of both ends. Then, a wide angular lens was used to take pictures at both extremes of the transect to achieve precise localization (usually +/- 1 m) in subsequent years (the same camera and lens should be used for this aim). Pictures can also be used to follow interannual changes in habitat structure ( Figure 1 ).

Habitat structure of sampled areas can be visually estimated with the help of reference charts (e.g. Prodon and Lebreton, 1981). We conducted an annual estimation of plant cover at the grass, shrub, subcanopy and canopy layers. To this aim, we sampled four virtual 5-m radius semicircles centred at 5, 35, 65 and 95 m from one end of the transect ( Figure 2 ) a few weeks after each exuviae sampling. If required for the study, more precise (and time-consuming) plant cover sampling using line transects or plots can be performed in sections of the exuviae transect.

Plots are used to describe the phenology of focal species in the study region (see under 4.6. Cicada emergence phenology) and determine sampling dates more precisely. They also allow densities of exuviae to be estimated and related with transect counts (see under 4.4. Estimation of absolute densities of exuviae). If cicada emergence phenology in the region is unknown, plot sampling has to be done before the start of transect monitoring. Four plots are located near a sample of transects, two at each end, following the design found in Figure 2 and Appendix B .

For this study we selected 15 transects evenly distributed among three vegetation age classes (one and two years after fire and unburned forests) and a total of 60 plots. Plots were 3 x 3 m squares in sparse vegetation (one and two years after fire) and 2 x 2 m in dense vegetation (>40 year-old forest). Exuviae were carefully searched by two people for four minutes per plot so that virtually all of them were collected. This procedure was repeated every 7 days across the summer, from before the onset of the emergence period until two weeks after the last exuvia was found. In our case, plots were sampled 13 times in 2014, from June 18 to September 11, finding four cicada species, Cicada orni being the dominant one (N= 315 C. orni of a total of 342 exuviae).

We tried to keep walking speed constant among observers, transects and years (around 12-15 minutes to walk 100 m without stops, sampling one side of the path). However, overall sampling speed was also affected by stops to collect exuviae and to measure field variables. In addition, sampling speed could be lower when plant cover was denser if more time was needed to revise potential moult supports. We used linear mixed models (LMM) with Gaussian error structure to analyse the effects of four explanatory variables (number of exuviae found, number of field variables measured, grass cover and shrub cover) on the sampling speed (m/min) of 283 transects. The ‘transect’, nested within ‘burned area’, was included as a random factor in order to control for possible spatial autocorrelation. Details on model selection can be found in Appendix C .

The best averaged model showed that sampling speed was affected by the number of exuviae found and the number of field variables measured, but not by grass and shrub cover ( Appendix C ). Both selected explanatory variables had a negative effect on speed, as predicted. However, the impact of finding a higher number of exuviae was more marked than the effect of increasing the number of variables measured on-site ( Figure 3 ), probably because the range of the latter was small.

A total of 24 transects (12 in burned areas and 12 in unburned areas) were used to assess detectability. They were sampled twice successively and at the same speed, with the second count progressing in the opposite direction to the first one. Exuviae were determined and counted as usual. We found 1,212 exuviae in the first walk of the 24 transects and 328 exuviae in the second. We analysed detectability for the two species with more than 20 records, Cicada orni and Tettigettalna argentata. We assumed that the proportion of undetected exuviae in successive counts (third, fourth…) would be the same as that calculated between the first and second counts. In C. orni, the second count produced 26.4% of the exuviae counted in the first walk (N=1166). For this species, two counts would already produce 93.0% of the total number of existing and visible exuviae. In T. argentata, the second count produced 46.9% of the exuviae counted in the first walk (N=32). For this species, two counts would produce 78.0% of the total number of existing and visible exuviae and four counts would be needed to reach 95.2% ( Table 2 ).

The detectability of C. orni exuviae appeared higher in burned sites than in unburned sites, with an estimated 77.5% vs. 69.2% of exuviae found in the first count, respectively ( Appendix D ). However, this difference was non-significant (t= -1.553, df= 22, p= 0.135). Exuviae numbers were insufficient to conduct the same burned vs. unburned comparison in T. argentata.

Exuviae fallen to the ground seemed less detectable than exuviae attached to vegetation. We used transects with at least 6 exuviae found at the second count (N=12), to compare the percentage of exuviae on the ground between the first and the second walk of the transect. Although this percentage was on average 13.6% in the first count and 21.2% second count (pointing towards fallen exuviae undetected in the first count because they were less visible), the difference between counts was non-significant (One-factor Anova, F=0.657, p=0.426).

We tested the effect of intense weather conditions on the persistence of exuviae on their supports, because fallen exuviae seemed harder to detect, especially under dense cover, and they were likely to degrade sooner. To this aim, we used 196 transect counts that met two requirements: having a minimum number of 15 exuviae, and having a recorded number of exuviae attached to vegetation or fallen from their support. We collected meteorological data from the closest weather station according to the Catalan Meteorological Service. Stations were located at distances ranging from 1.5 to 13.5 km from selected transects. We used half-hour data (48 periods per day) for the period encompassing from 1st July to the date of each count, for each year of sampling. The response variable analysed was the proportion of fallen exuviae out of the total number. As we hypothesized that exuviae were dragged down by the intensity of meteorological phenomena, we built a model using the maximum rainfall (mm) and the maximum wind speed (m/s) in a half-hour period, as explanatory variables, including also the quadratic term of the wind speed because of its non-linear relationship with the proportion of fallen exuviae. The study area was used as a random factor to account for the spatial correlation of transects.

Both variables were significant in models ( Appendix E ). The effect of rainfall, which had been intense in a small number of meteorological stations and emergence periods, was positive and clearer than that of the wind ( Figure 4 ). Our model predicted that, in absence of rain, less than 10% of exuviae would have fallen to the ground when the transects were covered, whereas for rainfall exceeding 20 mm in half an hour, half of them would have fallen. Wind had a positive quadratic effect, with no clear effects at low wind intensities and exuviae starting to fall at wind speed above 20 m/s. It is important to note that both wind and rain are usually mild in our study region during the emergence period.

The accumulated counts of exuviae in the study plots throughout the whole emergence period provided a precise estimation of the density of emerging final instar nymphs. In small plots, the proportion of missed exuviae was greatly reduced after a careful search. We found a variable density of exuviae in the three vegetation age classes and between two species, C. orni and T. argentata ( Table 3 ). For both species, the highest density was found one year after fire, reaching a mean density of 121 emerged C. orni nymphs in 100 m² across one season (N= 20 plots).

Linear models show that counts in the 15 selected transects were strongly related to both the density of exuviae in nearby plots until the date when the transect was sampled and the density of exuviae in the same plots for the whole emergence period. The equation (y=0.0994x + 0.3226) indicates that the transect counts accounted for 10.3% of the exuviae emerged in plots until the transect sampling date. Similarly, transect counts accounted for 9.9% of the exuviae emerged throughout that season (y=0.0985x + 0.0996) ( Figure 5 ).

The ability of Cicada-MET to record species occurrence can be compared to common methods used for the study of cicada distribution and ecology. Since adult cicadas use sound to communicate, bioacoustics has been a prominent field of study in Cicadidae. From 2014 onwards, we recorded the presence of singing males of cicada species, while walking along the transects. A total of 394 exuviae counts coupled with singing information of five species was available. Aural species determination on-site by trained observers was based on sound recordings and sonograms found in Gogala (2009). Most species were easily distinguished; only Tibicina garricola was recorded and photographed to confirm species identity. Each species could be either undetected, or detected only by exuviae, by sound or by both methods in each transect ( Figure 6 ). We used the paired-samples McNemar test to test for the difference in the frequency of detection exclusively by exuviae or by sound, for each species separately. C. orni was significantly more detected by exuviae (p<0.001). Cicadatra atra (p<0.001) and T. argentata (p<0.05) were significantly more detected by sound. Lyristes plebejus (p=0.32) and Tibicina garricola (p=0.22) showed non-significant differences.

The consecutive sampling of fixed plots across the emergence season allowed emergence phenology to be studied. To this aim, we used the 60 plots previously described, equally distributed among burned (one and two years after fire) and unburned pine forests. The temporal evolution of density across the 2014 season followed a unimodal pattern in the three vegetation types. However, the corresponding emergence peaks were separated by seven days each. The peak was around July 3 in one-year postfire forest, July 10 in two-year postfire forest, and July 17 in unburned forest ( Figure 7 ). The length of the period with cicada emergences also increased with increasing time since fire, from 5 weeks in one-year postfire forest to 10 weeks in unburned forest.

To determine if there were temporal differences in the number of captures between the three vegetation types, a repeated measures analysis of variance was used ( Appendix F ). Results showed that both time since fire and calendar date were relevant elements. Moreover, the interaction of both temporal variables indicated that the effect of time since fire on cicada numbers depended on the day on which the measurements were taken. In other words, temporal trends were different among vegetation categories.

Fixed transects repeated yearly allowed monitoring of population trends over longer periods of time. As an example, we used 21 transects distributed in three study areas of unburned Aleppo pine forest in L’Empordà region, sampled between six and ten years (2013-2022). The transects selected were separated by a maximum distance of 31 km, in similar altitudes (40-350 m a.s.l. and distance to the sea (15-26 km), so that climate effects were likely irrelevant. The abundance of emerging Cicada orni in the three study areas was different, especially when comparing more distant La Jonquera to the two closer areas of Colomers and Vilopriu ( Figure 8 ). However, interannual variation was very similar, the only exception being the period 2015-2017, even with a modest number of transects per study area.

We propose a conceptual diagram to synthesize the factors affecting the number of exuviae recorded during transect sampling and the direction of their influence ( Figure 9 ). Some of these factors have been analysed in the present study. The observer’s skill and speed of progression directly influence the number of exuviae that are found or left behind. Although no data was available to assess observers’ skill in our case, observer differences are a common bias in visual censuses (see, e.g., Bernard et al., 2013). The number of exuviae collected significantly influenced the time needed to cover a transect, because more stops were made. Conversely, walking speed necessarily affects exuviae numbers because higher speed means less time for properly searching vegetation and the ground, although an experimental approach would be needed to describe this relationship. Other drivers of exuviae detectability could have been included but their effects are poorly known. For example, some animals (ants, other insects and small vertebrates) might consume or carry exuviae but, to our knowledge, this has not been described and is unlikely to significantly alter counts (we have never observed this behaviour in our sites).

Weather has several implications on the results of the method. Soil temperature, when reaching around 18°C, is the proximate cue for last instar nymphs of studied species to emerge (Heath, 1968). In warmer pre-emergence periods or sites, the soil is heated since an earlier date, triggering earlier nymph emergence (Smith et al., 2006; Sato and Sato, 2015). Large interannual (and spatial) differences in emergence phenology could entail the sampling date to be sub-optimal. Sampling could then be too early, i.e. before the peak of emergence, detecting few exuviae; or too late, i.e. too many days elapsed since the emergence peak, so that many exuviae have fallen from their support, degraded or disappeared. If this is suspected, repeated plot sampling across the emergence period would be necessary. We have also shown that intense rain and wind increase the number of exuviae fallen from their support. Fallen exuviae can be dragged out of the sampled area or remain on-site in the ground and be found, although they are usually less apparent than when attached to vegetation. Detectability tends to decrease under these circumstances.

Habitat characteristics are important in exuviae counts. Plant cover did not influence observer speed, according to our results. It can, however, affect the detectability of exuviae since attached and fallen exuviae are less visible under high foliage density (pers. obs.). The shadow produced by plant cover also limits solar radiation and therefore soil temperature (Smith et al., 2006), with the consequences on emergence phenology discussed above. Interannual differences in host plant productivity or quality can also entail more or fewer nymphs being ready for emergence each year (Karban, 2022). Finally, spatial and temporal variations in habitat quality, predation, parasitism and competition influence the density of cicada species at each location and year.

In sum, the number of exuviae recorded in transects is not only determined by the number of emerging nymphs at each site and time. Multiple factors affecting the timing of emergence and the detectability of exuviae also influence the counts generated. This will be discussed in more detail in the next section.

Our results show that Cicada-MET is an effective and efficient method to count emerged cicadas in natural and semi-natural habitats. The amount of data obtained per time invested is high. In our study areas, we obtained an average of 24.9 exuviae (all species combined) per hour of transect (N= 350 counts, randomly selected, therefore including zero counts), while simultaneously measuring field variables related to location, habitat and moult support. Eight to twelve transect counts can be done per day by two observers working together, with adequate light conditions and avoiding the hottest hours of the day. Likewise, it is possible for two people to count 60 nearby plots in one day. The method needs minimal and non-expensive material and can therefore be used when research funds are scarce.

Sampling plots along the whole emergence period provide absolute densities of emerged cicadas. Our results show that they are very highly correlated to single transect counts, making it possible to estimate cicada densities from transects with measurable error. We estimated that transects capture about 10% of the total number of emerged nymphs in surveyed areas. Cicada-MET can therefore be used in ecological studies for which knowing absolute numbers of individuals or biomass of cicadas is necessary, for example, for understanding nutrient pools, trophic networks or biotic communities (Yang and Gratton, 2014).

The density of cicada nymphs along paths and tracks may be different to that of other parts of the habitat. This could result from females preferring sunnier locations in forest edges to oviposit (Yang, 2006), where trees grow more rapidly and attain greater root biomass (McDonald and Urban, 2004). Nymphs can also grow heavier and adult females be more fertile, in these areas, as indicated by slight differences in a 17-year life cycle species (Karban, 2014). On the other hand, the lack of vegetation on the track reduces oviposition rates in the transect or plot vicinity. Caution is therefore needed to infer mean densities for larger areas. A way to overcome this limitation can be to sample fixed plots at increasing distances from the path and compare the differences in exuviae abundance between them.

The method assumes that none or insignificant numbers of exuviae can remain from one year to the next (Moriyama and Numata, 2015, pers. obs.). If there are reasons to think that some proportion of exuviae can last on-site for one year, it would be necessary to remove them from transects and plots before the start of the emergence period.

The method should be adapted to study habitats holding the extremely high-density species that occur among cicadas (Karban, 1982). This can be achieved by reducing the area sampled, the number of variables measured, and/or by collecting only a sample of the exuviae found. If the study deals with periodical cicadas, in addition to these modifications, sampling has to be adapted to match the temporal and spatial occurrence of 13-year and 17-year species.

The standardised nature of Cicada-MET, a protocol that records precise location, area and transect speed, reduces the bias between sites that have been pointed out in former methods (Moriyama and Numata, 2015). The estimation of species detectability using repeated transects can allow interspecific comparisons in relative abundance of emerging cicadas. Additionally, the measurement of population densities with fixed plots makes such comparisons more robust. When possible, it is worth taking weather, observer differences and walking speed into account in models of exuviae abundance. Even if it is not, the results are highly consistent (see sections 3.2, 3.4, 3.5 and 3.7). To cope with emergence phenology variations in the longer term due to environmental change, plot counts across the emergence season could be repeated every 5 to 10 years.

Cicada-MET performs well when compared to auditory detection (presence/absence) of cicada species. Cicada orni, a loud species that can be heard from hundreds of meters away, was detected more frequently by exuviae than by hearing in our transects. To avoid the hottest hours of the day, transect sampling started at sunrise, when cicadas have not usually started to sing because of cooler temperatures. This can explain the better performance of exuviae detection, but is also a benefit for fieldwork in hot regions where sampling during the peak cicada song period can be a challenge. For other species, we obtained similar or more auditory detections. In Cicadatra atra, we recorded few individuals by exuviae counts because the pine forests surveyed are not a preferred habitat, whereas we could hear their song from nearby broadleaf trees. Nonetheless, the association between exuviae numbers and the abundance of recently emerged imagos is direct. In contrast, estimating animal densities from sound recording, albeit a promising field of study, needs expensive material and important data storage and processing capacity (Darras et al., 2019; Pérez-Granados and Traba, 2021).

Our results confirm the advance on peak emergence date in burned vegetation compared to unburned vegetation, which has also been found in riparian forests of North America (Smith et al., 2006). That study found earlier dates and shorter periods of emergence in burned forests related to mean soil temperature, being lower in sites with higher tree cover. Cicada-MET provides a fine-resolution description of phenology based on replicated plots. We specifically found a one-week delay on the peak emergence date between the first and second year after fire and between second-year vegetation and unburned forests. Our study therefore demonstrates that wildfire can have significant impacts on cicada emergence phenology, and highlights the importance of understanding how natural disturbances (Viljur et al., 2022) affect the timing and duration of population processes. Further research is needed to fully understand the ecological implications of advanced cicada emergence and how it may impact cicada mortality and reproduction as well as the derived interactions in ecosystems.

The similarity we found in interannual abundance trends between nearby study areas is remarkable ( Figure 8 ). As discussed above, this finding supports the consistency of the method. Transects are effective in capturing population numbers while biases and errors do not produce important random differences in temporal patterns between areas. The pattern found also leads us to assume that factors governing the quantity of emerging cicadas operate evenly at a landscape scale in the study region. Interannual climate variability seems to be the main factor here, because it affects plant productivity, disturbance regime, species interactions, cicada reproduction, and survival and length of the nymphal phase (Boulard, 1965; Williams and Simon, 1995), therefore modulating cicada numbers. On the other hand, poor nutrition (xylem fluid exclusively) and intrinsic biological traits also determine the length of the nymphal phase (Karban, 2022). Altogether, it appears that some interannual spatial synchronicity is possible in common species such as Cicada orni. In its extreme form, emergence synchronicity occurs in periodical cicadas, which include seven species in North America, one in India and one in Fiji (Simon et al., 2022).

C. orni was the most abundant species in all study areas and years. When considering fixed plots, C. orni mean densities were higher in recently burned (121 emerged nymphs/100 m²) than in unburned Aleppo pine forests (29 emerged nymphs/100 m²), across the emergence season. When considering the 466 transect counts, only 9.4% of them did not detect the species whereas the mean count was 11.7 exuviae/100 m². Although this figure corresponds to pine forest located along paths, which could hold denser nymph populations than more shadowed areas, it is within the range of former counts on transects located outside paths (Pons, 2015). The maximum recorded count in our transects was 191 exuviae/100 m², corresponding to an absolute density of 1,926 +/- 331 emerged nymphs/100 m² across the whole season, according to our plot-transect model. This is a considerable density, but still comparable to former estimations. Patterson et al. (1997) estimated densities between 712 and 1,117 of belowground C. orni final instar nymphs/100 m² in Italian pine forests (most of them emerging that year, according to their counts of emergence holes and exuviae). These numbers are, however, far from the densities found for periodical cicadas in North America, which can reach locally 57,900 emerged nymphs/100 m² (Karban, 2014).

Cicada-MET can obtain a large number of replicates and be used in a variety of ecological studies. This is facilitated because counting cicada exuviae can be easily combined with collecting diverse variables of interest. Measurements can be done on-site while sampling the transects or plots. They can be focused on exuviae ( Table 4A ), such as determining taxonomic identity, abundance, location, type and use of moulting support; or consider their fine-scale environment, including the lithology, topography, meteorology, and habitat composition, structure and function. Environmental data at a coarser scale can be obtained from Geographical Information Systems, remote sensing databases and UAV sensors ( Table 4B ). Due to the low mobility of nymphs, both belowground and once emerged heading to a moulting support (Oberdörster and Grant, 2006; Karban, 2014), the measurements in exuviae locations allow testing spatial hypotheses about oviposition and nymph occurrence sites (Bergh et al., 2022). The coupling of environmental and biological variables thus obtained increases the value of cicadas as soil and vegetation indicators.

Other research questions can benefit from the collection of exuviae in transects and plots of known location ( Table 4C ). Stored exuviae can last for long periods without degradation, more than one decade in absence of preservative (own collection at Univ. of Girona), and more than one century in museum collections with preservative (Cicada orni exuviae coded 75-8984, dated 1900 at the Natural Sciences Museum of Barcelona). Preserved exuviae can be used to determine sex ratios, to do biometric studies, or for genetic, isotopic or ecotoxicological studies. Their analyses can indicate population structure, dietary changes, root and soil health, etc. (Zheng et al., 2010; Lee et al., 2012; Bouwer et al., 2014).

Depending on the temporal and spatial arrangement of transects, Cicada-MET can be used at different scales going from intensive studies focused on habitats or study sites, to large-scale monitoring programs. The method is well suited for Before-After-Control-Impact studies (BACI) because it allows adequate replication in space and across time (Stewart-Oaten et al., 1986). It is also usable in ‘opportunistic’ studies after ecological disturbances in terrestrial ecosystems, including wildfires, wind blows, biological invasions, tree mortality due to insect explosions, or water stress. For this type of studies, the simplicity of the method facilitates an immediate start after the disturbance, and also a long-term approach (e.g. Figure 8 ), which is fundamental for population biology studies (Reinke et al., 2019).

Population densities derived from standardized counts of exuviae can be relevant for conservation purposes, as illustrated by dragonfly studies (Khelifa et al., 2018). Since cicadas are sensitive to changes in habitat quality and environmental condition, Cicada-MET data on species composition and abundance can produce informative ecological indicators ( Table 4 ). Changes in environmental conditions affect the timing and success of emergence (Smith et al., 2006) as cicadas are highly dependent on their host plants and the overall health of the soil. Emerged numbers over long periods can therefore serve as an indicator of changes in habitat quality and ecological integrity (Karban, 2014). Monitoring cicada populations and their emergence patterns therefore provides evidence to help act upon ecological conditions and long-term trends of terrestrial ecosystems.