We found that increased rainfall reduced exuvial numbers and impacted populations of Taiwanosemia hoppoensis from coastal regions, likely due to prolonged soil flooding harming the subterranean nymphs. As climate change advances, long-term monitoring is essential to track cicada populations across taxa and regions, given their important ecological roles.

Exuvial sampling and weather data

The study was conducted in the Chengxi windbreaks, dominated by Australian pines (Casuarina equisetifolia). We establishing four plots, spaced about 500 meters apart. In 2015 and 2016, we selected 12 same Australian pines per plot for exuvial sampling. In 2019 and 2022, we retained the original trees if they remained alive; if any had fallen or died, we replaced them with nearby Australian pines and increased the sample size to 15 trees per plot. Surveys were conducted between mid-April and mid-August of the selected year (45 surveys in total: 11 in 2015, 10 in 2016, 13 in 2019 and 11 in 2022). Each session occurred every 5 to 17 days and lasted about 4 hours. Exuviae were collected and counted by hand from the trees and within about a 1-meter radius on the ground. Although we occasionally found exuviae on the ground, their numbers were minimal.

We obtained daily air temperature, relative humidity, daily accumulated rainfall, and monthly rainy day counts from the Anping weather station (22°59'34" N, 120°09'07" E, about 9 km southeastern of the study area) for the period between August 2012 and July 2022. Although exuvial surveys began in 2015, weather data from 2012 onward provided a baseline for understanding pre-study weather conditions, particularly rainfall.

Emergence season scale

To assess the impact of weather conditions during emergence periods on exuvial numbers, we calculated mean air temperature (mean T_a), mean relative humidity (mean RH), total accumulated rainfall (TAR) and simple daily intensity index (SDII, TAR divided by the numbers of rainy days) between two exuvial survey periods (n = 45 for all). To account for irregular survey intervals, we standardized exuvial numbers, TAR, and SDII by dividing them by the number of days between surveys, yielding standardized numbers of exuviae (STD_NE), standardized TAR (STD_TAR), and standardized SDII (STD_SDII), respectively. The STD_NE were rounded to whole numbers for GLMM (generalized linear mixed model) analysis.

For GLMM analysis, we began by screening the variables through Spearman's rank correlation between STD_NE and each meteorological factor (mean T_a, mean RH, STD_TAR and STD_SDII). Relative humidity was excluded from the GLMM due to its low correlation with STD_NE (Spearman’s rho = -0.24, p = 0.12, n = 45). In the GLMM, we modeled STD_NE as a function of three fixed-effect variables (mean T_a, STD_TAR, STD_SDII) and one random-effect variable (year) (i.e., the full model). Given the count nature of our data, we used a Poisson family with log links and included observation-level random effects to address overdispersion (= 3.36) (Harrison, 2014). We excluded the interaction terms (mean T_a ´ STD_TAR, mean T_a ´ STD_SDII and STD_TAR ´ STD_SDII) to prevent multicollinearity, as these variables were significant correlated, which could lead to convergence issues and unreliable results (Izenman, 2013). Multicollinearity was further assessed using variance inflation factors (VIF) (Dormann et al., 2013), finding high multicollinearity between STD_TAR (VIF= 31.26) and STD_SDII (VIF= 31.84). We retained only STD_SDII in the final model due to its stronger correlation with STD_NE (Spearman's rho= -0.47, p = 0.001, n = 45) compared to STD_TAR (Spearman's rho= -0.38, p = 0.01, n = 45) (Figure S5). While only the final model results were reported in the main text, both models were presented in Table S1 for reference. The GLMM analysis was performed in R-4.1.1 using the lme4 package.

Interannual scale

We calculated total exuviae collected in 2015, 2016, 2019, and 2022, We also calculated one- and two-year values for mean T_a (MT_{a one-year}, MT_{a two-year}), one-year and two-year accumulated rainfall (A_TAR_one-year, A_TAR_two-year), and one-year and two-year simple daily intensity index (A_SDII_one-year, A_SDII_two-year) for the period from previous August (and the August of the two preceding years) to July across the years 2013/2014 to 2022 (e.g., A_TAR_one-year in 2015 represents the total accumulated rainfall from August 2014 to July 2015). We selected this timeframe because twilight cicada nymphs typically begin development in August (Chang et al., 2021), with a possible 2-3 year nymphal stage (or even shorter). Moreover, their emergence (early May to early August) coincides with Taiwan’s rainy season (Chang et al., 2021), when heavy rains impact both emerging adults and underground nymphs.

We also calculated one-year accumulated rainfall and SDII for the rainy (August – September, May – July) (A_TAR_{rainy season}, A_SDII_{rainy season}) and dry (October – April) (A_TAR_{dry season}, A_SDII_{dry season}) seasons during the exuviae collection periods (August 2014 – July 2015, August 2015 – July 2016, August 2018 – July 2019 and August 2021 – July 2022) as an example to assess the effects of dry and rainy seasonal rainfall on cicada emergence, despite their 2-3 year nymphal stage.

We calculated the percentage decrease in the annual total exuviae collected in 2016, 2019, and 2022 relative to the total counts in 2015. Changes in A_TAR_one-year/A_TAR_two-year and A_SDII_one-year/A_SDII_two-year from 2016 to 2022 were expressed as quotients, with each year's values divided by those from 2015, when our exuvial survey began. Similarly, changes in A_TAR_{rainy season}, A_TAR_{dry season}, A_SDII_{rainy season} and A_SDII_{dry season} for 2016, 2019, and 2022 were expressed as quotients, using the 2015 rainy season values as the baseline. Data were expressed as means ± 1 SE (standard error), with significance level (α) of 0.05.

References

Chang, Y. M., Jong, J. J. & Lin, S. F. (2021). Taiwan Twilight Cicada. Taijiang National Park Headquarters.

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., Gruber, B., Lafourcade, B. & Leitão, P. J. (2013). Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography, 36, 27-46.

Harrison, X. A. (2014). Using observation-level random effects to model overdispersion in count data in ecology and evolution. PeerJ, 2, e616.