Urban heat island effect on cicada densities in metropolitan Seoul

Weather data

Firstly, we obtained meteorological data of 2014 from the Korea Meteorological Administration (http://www.weather.go.kr/weather/observation/aws_table_popup.jsp, accessed 21 May 2015). Daily minimum and maximum temperatures of 38 automatic weather stations within and in the vicinity of metropolitan Seoul were compiled from June 1 to August 31, 2014 (92 days in total). Among 38 weather stations, eight were determined as surrounding stations based on their geographic locations outside of the administrative boundary of metropolitan Seoul. This period was chosen because important life history events of these two cicadas, such as mating and oviposition, occur in summer. Moreover, although cicadas overwinter in diapause (Moriyama & Numata, 2008) and UHI is most intense in winter (Kim & Baik, 2005), the severity of winter weather conditions was not reported to influence hatching success of cicadas (Moriyama & Numata, 2009).

Averages of minimum and maximum temperatures in 92 days of each weather station were considered representative of minimum and maximum temperature of that station, respectively. Minimum and maximum temperatures of urban areas were 21.02 ±  0.99 °C (Mean ± SD) and 29.03 ± 0.66 °C, respectively; those of neighboring areas were 20.16 ± 0.72 °C and 28.58 °C, respectively. Independent-samples t-tests revealed significantly higher minimum temperature of urban stations compared to that of neighboring stations (P = 0.029) but no significant difference in maximum temperature between urban and neighboring stations (P = 0.093) (Supplemental Information 4). Therefore, minimum temperature was employed as an estimate of ambient temperature. Minimum temperature is also used in other UHI research since the difference in maximum temperatures between urban and suburban areas usually indicates a cooling rather than heating effect (Tzavali et al., 2015).

Sampling localities

We randomly chose eight urban and four surrounding locations for sampling, each location corresponding to one weather station. Three localities in each location were further chosen as replicates. Localities were standardized as residential complexes where landscaping trees were usually found. Since all three localities within each location were examined thoroughly within one day in each sampling, the boundary of each locality was set within the geographic areas of about 10,000 m². Besides, the distance between radii of any two localities was greater than 100 m, to avoid overlapping cicada ranges (Karban, 1981). Each locality encompasses two to six buildings within a residential area. Detailed information of each locality is provided in Table 1.

We visited all localities in each location three times, with an interval between consecutive visits of at least 14 days. In the first visit, remaining exuviae of preceding broods were removed. The second and third visits respectively represented the first and second samplings of cicada emergence. In total, 36 localities of 12 locations were sampled twice from June 17 to August 3, 2015 (Fig. 1).

Cicada exuvia collection

Since there is 1:1 matching between abundances of exuviae and adult cicadas, an enumeration survey of exuviae is a good predictor of cicada population density at a site (Moriyama & Numata, 2015; Sato & Sato, 2015; Takakura & Yamazaki, 2007). At the last developmental stage, the nymphs of the final cicada instars emerge from underground to molt. They shed their cast skins, or exuviae, on a nearby tree or on an artificial structure, then the adults fly away. These leftover materials can persist for a long time in the environment, even under exposure to variable environmental conditions. Accordingly, we relied on an enumeration survey of cicada exuviae to estimate population densities of the two target cicada species. Exuviae were collected from trees, underneath leaves, and on tree branches. Species identification of each exuvia was based on Lee, Oh & Jang (2012). Sexes were also determined from exuviae based on the morphology of the abdomen (Sato & Sato, 2015). Since there were no significant differences between the number of males and females in two sampling periods regardless of species (P > 0.05, Supplemental Information 5), the total number of exuviae was employed for the calculation of cicada population density.

We computed two types of resource-weighted densities as measurement of number of individuals living or using a resource (Lewontin & Levins, 1989): (1) area-weighted density, which is the number of exuviae divided by the total geographic area in a locality; and (2) tree-weighted density, which is the number of exuviae divided by the number of trees in a locality, regardless of the existence of exuviae. Area-weighted density indicates the number of cicadas living in a geographic unit, while tree-weighted density indicates the number of cicadas per available resource unit. Cicadas largely depend on trees during their lifetimes; they spend their nymphal stage underground, sucking fluid nutrition from tree roots; after emergence, imagoes continue relying on trees for food or as substrates for mating and oviposition. Consequently, trees can be considered as a resource for cicada population to exploit.

Environmental variables

To determine environmental factors that may predict cicada density, we used the values of wetness, greenness, and imperviousness inferred from satellite data, as well as minimum temperature values previously collected from weather data. Wetness can be correlated to soil and vegetation moisture (Lobser & Cohen, 2007), while greenness represents the measure of “photosynthetically active vegetation”, and can, therefore, serve as a proxy for primary productivity (Lobser & Cohen, 2007). Imperviousness consists of building footprints as well as pavement and asphalt.

Greenness, wetness, and imperviousness were derived from corrected reflectance Landsat 8 OLI/TIRS layers 2–7 using tasseled cap transformations (Supplemental Information 6; Baig et al., 2014; Grant & Carter, 2011). The formula for imperviousness was modified from the normalized difference vegetation index as an index for urbanization. Landsat layers were provided by the United States Geological Survey (USGS) and represent reflectance imagery from late May of 2014–2016.

TOA Reflectance is given by: ${\displaystyle \rho {\lambda }^{\prime }=M_{\rho }{Q}_{\mathrm{cal}}+A_{\rho }}$ in which ρλ′ is Landsat band calculated for surface reflectance, M_ρ reflectance multi-band x (from metadata), A_ρ reflectance add band x (from metadata), and Q_cal original Landsat band.

TOA Corrected Reflectance is calculated as: ${\displaystyle p\lambda =\rho {\lambda }^{\prime }∕sin\theta \mathrm{SE}}$ where ρλ corresponds to reflectance band corrected for sun elevation, and θSE sun elevation (from metadata).

Imperviousness is measured by: ${\displaystyle \mathrm{Imperviousness}\left(\mathrm{index}\right)=1-\frac{\mathrm{Band}5-\mathrm{Band 4}}{\mathrm{Band 5}+\mathrm{Band 4}}.}$

The greenness, wetness, and imperviousness values were then extracted to survey points to use in linear regression models.

Statistical analysis

A non-parametric Kruskall–Wallis one-way ANOVA test was used to compare species composition across sampling localities.

Multiple linear regression was carried out to determine significant factors for resource-weighted densities of two cicada species. The dependent variable was either area- or tree-weighted density of each species. Independent variables consisted of sampling period, minimum temperature from June to August 2014, age of residential complex, wetness, greenness, and imperviousness. Independence of residuals was examined by a Durbin-Watson test. A Durbin-Watson statistic ranging from 0 to 4 indicates no correlation between residuals. Multicollinearity was inspected by correlation coefficient and variance inflation factor (VIF). None of the correlation coefficients were greater than 0.7 (Supplemental Information 7), and all VIF values were less than 3, which confirmed no multicollinearity. Normality of the residuals was assessed visually by Q–Q plot of studentized residuals. Resource-weighed densities of both species required transformation to meet normal distribution assumption of residuals: for C. atrata, area-weighted density: log(x + 1), tree-weighted density: log(1, 000x + 1); for H. fuscata, area-weighted density: log(x + 1), tree-weighted density: log(10, 000x + 1), in which x  is the density of each species. Homogeneity of variance of residuals was examined using diagnostic plots of predicted values versus standardized residuals. Significant outliers were investigated using Leverage values and Cook’s distance. All Leverage values and Cook’s distances were less than 0.3, which suggested an absence of significant outliers. Overall, all assumptions of linear regression models were met. All statistical analyses were performed in SPSS 22 (IBM Corp; Armonk, NY, USA).

Resource-weighted densities of C. atrata

Area-weighted density was 119.71 exuviae/km²/day, and tree-weighted density was 3.1 × 10⁻³ exuviae/tree/day in the first sampling period. The resource-weighted densities in the second sampling was triple those in the first sampling. Overall, densities of C. atrata increased with the increase in minimum temperature.

Linear regression models significantly predicted resource-weighted densities of C. atrata (area-weighted density F_6,65 = 9.76, P < 0.001, adjusted R² = 0.42, tree-weighted density F_6,65 = 9.19, P < 0.001, adjusted R² = 0.41). Sampling period and minimum temperature were significantly positively correlated to both area- and tree-weighted densities of C. atrata (P < 0.05, Table 2, Fig. 3), whereas age of residential complex and other environmental factors were not significant (P > 0.05, Table 2).

Resource-weighted densities of H. fuscata

In the first sampling, resource-weighted densities were respectively 132.18 exuviae/km² /day for area- and 4.72 × 10⁻³ exuviae/tree/day for tree-weighted density. Compared to those, the resource-weighted densities in the second sampling were nine times higher. Analogous to C. atrata, the overall pattern represents a parallel increase between H. fuscata densities and elevated minimum temperature.

Linear regression models were significant for H. fuscata densities (area-weighted density F_6,65 = 18.93, P < 0.001, adjusted R² = 0.60, tree-weighted density F_6,65 = 19.97, P < 0.001, adjusted R² = 0.62). Again, sampling period and minimum temperature were significant variables for resource-weighted densities (P < 0.001, Table 3, Fig. 4). In addition, age of residential complex did not exhibit a significant correlation to C. atrata densities, but to H. fuscata densities (P < 0.001, Table 3, Fig. 5). Besides minimum temperature, no other environmental factors were found significant (P > 0.05, Table 3).