Not all cicadas increase thermal tolerance in response to a temperature gradient in metropolitan Seoul

Rapid anthropogenic alterations caused by urbanization are increasing temperatures in urban cores, a phenomenon known as the urban heat island (UHI) effect. Two cicada species, Cryptotympana atrata and Hyalessa fuscata (Hemiptera: Cicadidae), are abundant in metropolitan Seoul where their population densities correlate strongly with UHI intensities. Such a positive correlation between cicada density and UHI intensity may be possible if cicada abundance is linked to a certain thermal tolerance. We tested this hypothesis by investigating variation in morphology and thermal responses of two cicada species along a thermal gradient in Seoul and surrounding areas. The morphological responses were measured by metrics such as length, thorax width and depth, and mass. The thermal responses were measured in terms of minimum flight temperature, maximum voluntary temperature and heat torpor temperature. First, we observed a species-specific variation in thermal responses, in which C. atrata displayed a higher thermal threshold for maximum voluntary and heat torpor temperatures than H. fuscata. Second, a positive association between temperature conditions and body sizes were displayed in females H. fuscata, but not in either conspecific males or C. atrata individuals. Third, C. atrata exhibited similar thermal responses regardless of habitat temperature, while H. fuscata in warmer areas showed an increase in heat tolerance. In addition, H. fuscata individuals with bigger thorax sizes were more heat-tolerant than those with smaller thoraxes. Overall, our research is the first to detect a variation in thermal responses and body size of H. fuscata individuals at a local scale. More investigations would be needed to better understand the adaptation mechanisms of insects linked to UHI effects.


Methodology
Sample collection. We sampled populations of C. atrata and H. fuscata from July 15 to August 5, 2016, in metropolitan Seoul and the vicinity in the Republic of Korea. Metropolitan Seoul covers more than 600 km 2 with diverse landscape features, generating a mosaic of heterogeneous UHI intensities for cicadas. Approximately 10 million people live within city limits, with another 10 million in the surrounding suburban areas.
The sampling method and selection of study areas followed 12 sampling areas by Nguyen, et al. 37 . From there, three areas were excluded due to low sampling densities and difficulty in collecting cicadas. In order to verify our sampling design and its suitability to our research question, we compared those nine areas in terms of the abiotic factors related to urbanization such as greenness, wetness, imperviousness and isothermal. Greenness is a measure of photosynthetically active vegetation, wetness represents soil and vegetation moisture 38 , and imperviousness measures building footprints, pavements and asphalt. Isothermal is mean diurnal range divided by temperature annual range, extracted from worldclim 1.4 39 . A one-way analysis of variance in which greenness, wetness, imperviousness, and isothermal were response variables was conducted. As expected, isothermal was the only variable that was significantly different among the nine areas (Supplementary Material 1). We concluded that temperature was the critical factor that differentiated among those nine area, thus verifying our sampling scheme (Fig. 1).
Weather data were compiled from Korea Meteorological Administration records. We calculated average maximum ambient temperature (T a Max) of each area during summer periods from June 1, 2010, to August 31, 2015. Given that mean, minimum, and maximum ambient temperatures are employed in UHI studies 40 , we relied on maximum ambient temperature, as it better represents the actual high ambient temperature that the cicadas encounter in their environment than the mean or minimum ambient temperature. Nine areas were randomly sampled and each area was sampled twice. Cicadas were collected from 8:00 a.m. to 12:00 p.m. at residential complexes in each area, and were subjected to thermal-response experiments within the day of capture. Information regarding sampling areas and the number of each species collected at each area are provided in Supplementary Material 2.
Measurement of thermal responses. Minimum flight temperature (MFT), maximum voluntary temperature (MVT), and heat torpor temperature (HTT) are conservative measures of thermal adaptation of cicadas to a habitat 4,41,42 (see 43 for a summary of thermal responses in 118 taxa of North American cicadas). MFT represents the lowest body temperature with fully coordinated activity. MVT is the upper thermoregulatory point at which thermoregulation takes precedence over other behaviors. HTT is the upper limit beyond which cicadas sink into a state of torpor. MVT increases as a habitat becomes warmer in some cases, whereas HTT is strictly related to the thermal condition of a habitat 4,[41][42][43][44] . The MFT-to-HTT range indicates the fully active thermal breadth (T b Range) of each species in relation to a certain habitat condition.
Thermal responses of each individual were assessed in a single assay with no rest time between treatments under laboratory conditions. First, each cicada was cooled to a torpid state by keeping it at −20 °C. We checked each individual's T b every three minutes to ensure they did not freeze, as this might affect HTT 42 . As soon as the individual was torpid, we assessed MFT by dropping the insect from a height of 2.5 m. If the insect could not perform the expected behavior, it was allowed to gradually warm up at ambient temperature for one minute before being reexamined. After that, the insect was placed under a heat lamp to obtain MVT. The T b at which the individual moved away from the heat source and started to seek shade was determined to be its MVT. The cicada was continuously heated under the heat lamp until no movement was observed, at which point its T b indicated HTT. The procedure was not lethal, as individuals could recover to normal active conditions after several minutes. We acknowledged that this assay may have induced stress to some extent in cicadas. However, such measurement of thermal responses have been performed over several decades 43 , and further investigation is necessary to assess tentative influences of this assay on the thermal performance of cicadas.
The evaluation of T b in cicadas is commonly conducted inside the mesothorax 43 . Here, we chose to assess T b from both the mesothorax and the pronotum. Indeed, pronotom may represent an evaporative cooling site 45 . Evaporative cooling is a key physiological thermoregulation mechanism in cicadas, as it provides a major cooling effect by dissipating excess heat [45][46][47] and facilitates the cicadas' endurance of high ambient temperatures 4 . The temperature of the pronotum therefore may signal a critical thermal threshold for the individual to regulate T b within its functionally active range and adopt necessary thermoregulation strategies to prevent excessive increase in T b .
All temperature measurements were performed using a digital thermometer with a k-type thermocouple (Omega; model #: 450-AKT; Norwalk, Connecticut, USA) sensitive to ± 0.25 °C. The total live body mass of each individual was determined using an Adventurer Pro Analytical (Ohaus; New York, USA) scale sensitive to ± 0.0001 g. We also measured body length, mesothorax width, and mesothorax depth using Digital Calipers (Insize Co., Ltd.; Georgia, USA) sensitive to ± 0.02 mm.
Statistical analysis. Comparative thermal responses of C. atratra and H. fuscata. A first constrained multivariate analysis, i.e., redundancy analysis 48 (RDA), was performed to compare thermal responses of MFT, MVT, HTT and T b Range between C. atrata and H. fuscata. Analyses were conducted separately for temperature measurements from the pronotum and mesonotum. C. atrata and H. fuscata responded differently to heat experiments. Therefore, we conducted hereafter analyses separately for each species. Variation in morphological characters of each cicada species. First, we applied an RDA to evaluate the responses of the morphological characteristics measured by total mass (mass), body length (length), mesothorax width (width), and mesothorax depth (depth) to sex and T a Max. Second, intersexual morphological differences were assessed by performing t-tests for normally distributed data of mass and Kruskal-Wallis tests for non-normally distributed data of length, width and depth. Furthermore, we also conducted linear regressions to examine the effect of temperature on each of the morphological characteristics, separating analysis for males and females. We examined the assumption of homogeneity of variance of residuals of each linear regression model via a diagnostic plot of predicted values versus standardized residuals.
Variation in thermal responses of each cicada species. A third RDA was used to assess the thermal responses of MFT, MVT, HTT, and T b Range by sex, T a Max and width. We performed independent RDAs for temperature measurement from the pronotum and mesonotum. Again, we compared intersexual differences in terms of thermal responses, employing t-tests for normally distributed data (MVT, HTT and TB) and Kruskal-Wallis tests for non-normally distributed data (MFT). Linear regressions tested how both male and female thermal responses, measured from pronotum and mesonotum, were influenced by T a Max and width. Diagnostic plots of predicted values versus standardized residuals were visualized to assess the assumption of homoscedasticity of residuals for each linear regression model. RDA results showed no significant thermal responses of C. atrata in the measurement from the pronotum or mesonotum, and it was exempted from follow-up intersexual variation tests.
Intersexual variation in thermal responses of H. fuscata. Finally, an RDA was performed to quantify intersexual differences in thermal responses according to T a Max and width. The Vegan package 49 on R Studio (Version 1.0.143) was used for all multivariate analyses, and the statistical significance of the entire model for each variable (marginal test) was evaluated using Monte-Carlo permutation tests (n = 999). Linear regressions were performed with SPSS 22 (IBM Corp.; New York, USA). Results are presented as the mean ± standard deviation.
Ethics declaration. Cicadas are common species in Republic of Korea. Neither C. atrata nor H. fuscata was listed as protected or endangered species in the "List of wildlife species prohibited for collection" issued by the Ministry of Environment, Republic of Korea and in the IUCN Red List. Therefore, no field permit was required for this study.

Comparative thermal responses of C. atratra and H. fuscata. Thermal responses measured
at the pronotum (n = 158) showed that C. atrata became fully coordinated at an MFT of 24 In general, C. atrata tolerated heat better than H. fuscata in terms of MVT and HTT. The RDA models explained 11.77% and 13.17% of the total variation in temperature measurement from the pronotum and the mesonotum, respectively (Table 1). Both showed that thermal responses were significantly different between C. atrata and H. fuscata. Indeed, for both the pronotum and mesonotum temperature measurements, the species factor was the most significant (p = 0.001), representing the first axis for both ordination diagrams ( Fig. 2A,B) and accounting for 80.6% and 76.33% of the inertia, respectively. Among four thermal responses, only MFT was greater for H. fuscata, whereas the others tended to be higher for C. atrata ( Fig. 2A,B).
Variation in morphological characters of each cicada species. Morphological  www.nature.com/scientificreports www.nature.com/scientificreports/ depth of 14.61 ± 0.63 cm, and a mesothorax width of 16.59 ± 0.75 cm. The result of the RDA showed that the effect of T a Max and sex did not significantly explain the morphological measurement matrix (Table 2).
H. fuscata (n = 258) was characterized by a total live mass of 1.60 ± 0.25 g, with a body length of 32.54 ± 1.80 cm, a mesothorax depth of 12.06 ± 1.05 cm, and a mesothorax width of 13.84 ± 1.24 cm. The RDA model explained 11.5% of the recorded variation, in which sex was the principal factor (p = 0.001), represented along the first axis and accounting for 88.72% of the variance (Table 2, Fig. 3). T a Max was also significant (p = 0.002) and was represented along the second axis, accounting for 11.28% of the variance. Females were significantly heavier than males, but shorter in total body length (Supplementary Material 3). H. fuscata individuals exhibited a decrease in mass (Fig. 4A) but an increase in mesothorax size relative to T a Max (Fig. 4B,C). Linear regressions showed that T a Max influenced male and female morphological characteristics differently (Table 3); T a Max was significant for only the mass of males but not for other characteristics; on the contrary, T a Max was significant for almost all characteristics of females, except length.
Variation in thermal responses of each cicada species. Although C. atrata individuals exhibited some changes in their thermal responses, the species responded similarly to heat regardless of habitat conditions. The results of the RDA showed that both sex and T a Max had no significant effect on thermal responses obtained from the pronotum (p > 0.05, Table 4). For the mesonotum, sex was the only significant factor explaining the thermal responses (p = 0.029).
The RDA model showed that thermal responses from the pronotum were mostly influenced by T a Max (p = 0.001), followed by sex (p = 0.011) and width (p = 0.011) ( Table 4). The model explained 11.76% of the total inertia, of which 87.32% was explained by axis 1 and 12.57% by axis 2 (Fig. 5A). Regarding the mesonotum, the thermal responses were mostly significantly influenced by sex (p = 0.001) and width (p = 0.001), followed by T a Max (p = 0.04) (Table 4, Fig. 5B). The thermal responses measured at the mesonotum were driven primarily by width rather than other factors (Supplementary Material 4).

Intersexual variation in thermal responses of H. fuscata.
For measurement from the pronotum, males and females had similar MFT and MVT. However, females tolerated heat significantly better (Fig. 6A) and had wider thermal ranges than males (Fig. 6B). Both sexes significantly increased their HTT and T b Range as T a Max increased (Supplementary Material 5), but no relationships were found between T a Max with either MFT or MVT (Table 5). Width displayed a significant negative effect on MFT of both sexes (Supplementary Material 6), while this factor was significantly positively correlated with HTT and T b Range of both males and females (Fig. 7).

Discussion
In sum, our analyses suggest a local adaptation of thermal responses and thorax sizes in H. fuscata populations distributed along a thermal gradient in metropolitan Seoul, but not for the other cicada species. Although the results of our redundancy analysis showed high values of unconstrained variance (>85% in all of them), significant differences in thermal responses between two cicadas and within H. fuscata were determined. Specifically, thermal responses of C. atrata measured by MFT, MVT, HTT and T b Range indicated a better tolerance to heat stimuli than H. fuscata. Furthermore, no relationship was observed between the ambient temperature of the habitat and either morphology or thermal physiology of C. atrata. On the contrary, H. fuscata from habitats with higher ambient temperature had substantially enlarged thoraxes, endured heat better, and held wider fully active thermal ranges.
Although ambient temperature was significant for the pronotum's HTTs and TBs of both sexes (Table 5), it was marginally significant for the mesonotum's TB of males and HTT of females (Supplementary Material 3). To justify the effect of ambient temperature on thermal responses of H. fuscata, we conducted the RDA and linear regression again, using Mass instead of Width as a covariate (results not shown). RDA result shows that T a Max was significant for thermal responses of H. fuscata regardless of pronotum or mesonotum temperature. Furthermore, the result of linear regression analysis displays a consistent result between pronotum and mesonotum temperatures. We therefore conclude that overall H. fuscata increased their thermal tolerance in accordance to the increase in ambient temperature.
The increase in heat tolerance of H. fuscata resembles other research on urban-adapted insects, which indicates a close association between thermal tolerance and localized thermal clines 29,30,33 . Research on thermal responses of cicadas across a wide geographic range provides evidence that cicadas are more tolerant of warmer environments [41][42][43] . Cicadas of 38 species inhabiting Mediterranean habitats display an elevated HTT in accordance with the local thermal characteristics, regardless of taxonomic position or the diversity of particular plant species 50 . Our results further imply a localized thermal acclimatization of H. fuscata.
Here, we observed interspecific differences in thermal responses between C. atrata and H. fuscata, in which the warmer the habitat, the greater thermal responses C. atrata exhibited compared to the other species. This disparity can be explained by variation in species geographical origins: C. atrata originated in subtropical regions, whereas H. fuscata originated in tundra regions 35 . Additionally, segregation in microhabitat niches may contribute to how each species utilizes its habitat for thermoregulation. C. atrata perches mainly on top of the canopy, where it is exposed to solar radiation, whereas H. fuscata is found mainly on tree trunks in shaded environments 51 . As a result, adaptation to individual thermal regimes has led to variation in thermal responses between these species, a pattern that is well-discerned in other cicadas inhabiting tropical habitats 52 .    www.nature.com/scientificreports www.nature.com/scientificreports/ In line with other studies on thermal responses of cicadas, we found HTT depends strongly on the maximum environmental thermal regime 43,44,53 . This positive relationship was observed in H. fuscata, but not in C. atrata, regardless of habitat conditions. Such difference may be partially due to the origins of the two species. C. atrata is, therefore, more prone to experiencing higher thermal regimes and is adapted to high thermal conditions in metropolitan Seoul, thus exhibiting no difference in thermal tolerance across heterogeneous ambient temperatures. Environmental constraints applied to populations of the acorn ant, Temnothorax curvispinosus, are greater in lower latitudes, causing a reduction in evolutionary thermal responses relative to populations at higher latitudes 29 .
In contrast to C. atrata, distributions of H. fuscata at higher latitudes expose this species to colder environments, and the warmer conditions of metropolitan Seoul may induce a thermal acclimatization to warmer temperatures. Not only does urban warming seem to select for thermophilic species, but it also extends thermal tolerance ranges by elevating heat tolerance 30 . Here, as a function of ambient temperature, the thermal range of H. fuscata was extended toward warmer habitats. Better heat tolerance and wider thermal active ranges promote colonization of microhabitat niches generated by urbanization.
Besides thermal responses, our study suggest contrast relationships between thermal conditions and morphological characters of H. fuscata females, while thorax sizes increased as ambient temperature increased, total mass  www.nature.com/scientificreports www.nature.com/scientificreports/ decreased. According to Bergman size clines, warmer environments usually trigger the growth rate of ectotherms. As a results, those from warmer environments tend to be bigger than the ones from colder environments 54 . The decrease of thorax sizes of females H. fuscata from warm to cooler habitats in this study shows support to Bergmann size clines. Besides, females cicadas developed at cooler environments were heavier than those at warmer environments. The flies examined by Crill et al. 55 also exhibited similar changes in body dry mass to our  Table 5. Linear regression analysis to assess the effect of T a Max and Width on thermal responses of males and females H. fuscata obtained from the pronotum. www.nature.com/scientificreports www.nature.com/scientificreports/ study. Accordingly, only female flies developed at higher temperature were heavier in terms of dry mass whereas male flies were insensitive to their developmental temperature. However, the underlying mechanism is unclear. In our study, we observe such an inverse effect of temperature on thorax size and mass of females H. fuscata. The evolutionary explanation for this phenomenon remains elusive, thus we aim to elucidate it in further study. Besides, more investigations are needed to better understand this phenomenon and to clarify if these first results are linked to an exceptional population of small individuals in the coolest areas or if this could be a general trend.
Insect metabolic activity is accelerated under exposure to higher temperatures, which triggers an increase in body size 56 . Furthermore, warming climates have been proposed as causes of escalating growth rates of cicada nymphs underground, resulting in larger body sizes within a fixed development period 22 . Females of the scale insect Melanasimpis tenebricosa inhabiting warmer tree canopies exhibit larger body size than those living in cooler tree canopies 57 . An increase of 2.5 °C in rearing temperature promotes the development of larger females in the southern green stink bug, Nezara viridula 58 . There is also evidence that larger body size may be driven by local adaptation to warmer habitats. In periodical cicadas, for instance, the more southerly diverged M. tredecim experiences higher temperatures as adults possess larger bodies 22 . Larger bodies provide multiple ecological and evolutionary advantages by contributing positively to survival, fecundity, and mating success 59 . As a result, a larger size could result from selection for greater fecundity in female cicadas 23 .
The relationship between body size and thermal tolerance has been explored in various terrestrial insects. Although there are cases where superior heat tolerance is found in smaller 60,61 or intermediate-sized insects 62 , larger species tend to tolerate higher temperatures better than smaller species do 57,58,62,63 . In our study, larger H. fuscata individuals had higher heat tolerance. Consequently, better thermal tolerance capacity is suggested to provide higher fitness to female cicadas 3 . Although body size of male H. fuscata was not correlated with ambient temperature, it was found to be associated with acoustic properties of cicada songs [64][65][66] , which act as premating signals to attract conspecific mates. Furthermore, songs of warmer males are higher in intensity and are able to travel further in the air 67 , increasing the transmission of their mating signals to prospective mates.

conclusion
Our thermal tolerance experiments indicate a local adaptation of thermal responses and thorax sizes along a thermal gradient of H. fuscata in metropolitan Seoul. Whether such variation in heat tolerance is caused by phenotypic plasticity or evolutionary adaptation to environmental conditions is unclear. However, acclimatization to anthropogenic perturbation as a consequence of urbanization may be partially responsible. This is the first study to notice variation in thermal tolerance of cicadas at the local urban scale. Our research highlights the importance of taking localized thermal regimes into consideration when examining species-specific responses to escalating urban warming caused by urbanization.