Predictability of temporal variation in climate and the evolution of seasonal polyphenism in tropical butterfly communities

Phenotypic plasticity in heterogeneous environments can provide tight environment‐phenotype matching. However, the prerequisite is a reliable environmental cue(s) that enables organisms to use current environmental information to induce the development of a phenotype with high fitness in a forthcoming environment. Here, we quantify predictability in the timing of precipitation and temperature change to examine how this is associated with seasonal polyphenism in tropical Mycalesina butterflies. Seasonal precipitation in the tropics typically results in distinct selective environments, the wet and dry seasons, and changes in temperature can be a major environmental cue. We sampled communities of Mycalesina butterflies from two seasonal locations and one aseasonal location. Quantifying environmental predictability using wavelet analysis and Colwell's indices confirmed a strong periodicity of precipitation over a 12‐month period at both seasonal locations compared to the aseasonal one. However, temperature seasonality and periodicity differed between the two seasonal locations. We further show that: (a) most females from both seasonal locations synchronize their reproduction with the seasons by breeding in the wet season but arresting reproduction in the dry season. In contrast, all species breed throughout the year in the aseasonal location and (b) species from the seasonal locations, but not those from the aseasonal location, exhibited polyphenism in wing pattern traits (eyespot size). We conclude that seasonal precipitation and its predictability are primary factors shaping the evolution of polyphenism in Mycalesina butterflies, and populations or species secondarily evolve local adaptations for cue use that depend on the local variation in the environment.


| INTRODUC TI ON
Phenotypic plasticity is the capacity of a genotype to produce different phenotypes in different environments. In heterogeneous environments, adaptive phenotypic plasticity can maintain high fitness by providing a tighter environment-phenotype matching than a canalized response can (Nylin & Gotthard, 1998;Reed et al., 2010;Via et al., 1995). However, phenotypic plasticity is only adaptive when some environmental cues(s) can reliably signal the forthcoming selective environment, thereby allowing an organism to use current environmental information to develop a phenotype that will have high fitness in the future environment (Bonamour et al., 2019;Moran, 1992;Reed et al., 2010;Tufto, 2015). In highly stochastic environments, or in stable environments where maintenance of plasticity is costly, the selection is expected to favour a reduced or complete loss of plasticity (DeWitt et al., 1998;Leung et al., 2020;Reed et al., 2010).
The adaptive nature of phenotypic plasticity is evident in seasonal environments, which fluctuate temporally in a cyclical manner and provide cues (e.g. temperature or photoperiod) that can reliably signal the forthcoming season (e.g. Bradshaw & Holzapfel, 2007;Brakefield & Reitsma, 1991). Moreover, animals in seasonal environments are typically confronted with strong alternating selection pressure for reproduction during a favourable period and survival during an unfavourable period (Kivelä et al., 2013;Tauber et al., 1986;Varpe, 2017). In such discrete environments, seasonal polyphenism, which is an extreme case of phenotypic plasticity (Nijhout, 1999;Shapiro, 1976), can provide a tight environment-phenotype matching by expressing distinct developmental pathways, and therefore distinct alternative forms in different seasons, depending on the cues experienced during a critical period of juvenile development (Beldade et al., 2011;Nijhout, 1999Nijhout, , 2003Smith-Gill, 1983). It has become increasingly apparent that the phenotype involved in such polyphenic responses extends from morphology to physiological, behavioural and life-history traits (Simpson et al., 2011).
The occurrence of polyphenic thresholds is one of the widespread mechanisms in which a continuous environmental input is translated into a switch-like response (Nijhout, 1999(Nijhout, , 2003Smith-Gill, 1983). For example, in dung beetles, only males that exceed the threshold body size develop long horns on their heads, whereas males below this threshold do not or develop rudimentary horns (Moczek, 1998;Moczek & Emlen, 1999). Alternatively, the underlying reaction norm can be more linear but the discrete nature of the environmental variation results in a switch-like response in the field with few, if any, intermediate phenotypes (e.g. Brakefield & Reitsma, 1991). When environmental cues vary across the geographic distribution of a species, populations may show local adaptations in polyphenism by evolving a new polyphenic threshold (e.g. Bradshaw & Holzapfel, 2001;Lindestad et al., 2019;Moczek & Nijhout, 2003). The evolution of polyphenic responses (e.g. diapause induction) of some insects to regional differences in day length and temperature in temperate regions provide examples of this phenomenon (e.g. Bradshaw & Holzapfel, 2001;Lindestad et al., 2019).
Here, we investigate how marked variation in precipitation and temperature and their predictability shapes the evolution of seasonal polyphenism in three communities of tropical Mycalesina (Nymphalidae: Satyrinae) butterflies and populations of a distantly related satyrine butterfly (Melanitis leda). Seasonal precipitation in the tropics results in distinct environments (i.e. wet-and dry seasons), which impose contrasting selection pressures across seasons and is the major driver of morphological and life-history evolution in Mycalesina butterflies (Brakefield & Reitsma, 1991, Roskam & Brakefield 1999, Brakefield et al., 2007. As an adaptation to seasonality, many Mycalesina species exhibit seasonal polyphenism with distinct Wet and Dry Season Forms (WSF and DSF, respectively) in the respective seasons (Brakefield & Larsen, 1984;Brakefield & Reitsma, 1991; Figure 1). WSF adults have conspicuous "eyespots" on the wing margins, which can function in deflecting predator attacks away from vital parts of the body, whereas DSFs lack eyespots and have a cryptic wing pattern, which provides better background matching when perching on the ground amidst the dry leaf litter (Brakefield & Larsen, 1984;Lyytinen et al., 2004;Prudic et al., 2015; Figure 1). WSF butterflies reproduce without diapause, whereas DSFs generally undergo reproductive diapause (Braby, 1995;Halali et al., 2020). Furthermore, the seasonal forms differ in a suite of traits such as body size, investment in reproduction, abdominal fat content (e.g. van Bergen et al., 2017), longevity (e.g. Pijpe et al., 2008) and behavioural traits (e.g. Prudic et al., 2011).
Studies, mostly confined to Mycalesina species from Malawi in East Africa, suggest that these butterflies use temperature as a cue (Brakefield & Reitsma, 1991;Kooi & Brakefield, 1999;Roskam & Brakefield, 1996;Windig et al., 1994) (Brakefield & Reitsma, 1991;Windig et al., 1994;Figure 2). Such marked seasonal fluctuation in temperature could, therefore, be used as a cue to predict the forthcoming season (Brakefield & Reitsma, 1991;Windig et al., 1994). This field observation is corroborated by common garden experiments, which show that rearing larvae at <21°C results in the DSF and >25°C results in the WSF (e.g. van Bergen et al., 2017;Brakefield & Reitsma, 1991;Kooi & Brakefield, 1999). However, many regions in the tropics where Mycalesina butterflies occur show a narrow variation in the mean monthly temperature over the year, and also a lack of field studies, except in East Africa, has hindered investigating whether temperature still acts as a major cue (see Roskam & Brakefield, 1999).
We sampled Mycalesina butterflies in two highly seasonal locations, Zomba in Malawi and Sanguem in India, and one much less seasonal location (henceforth aseasonal), Genting in Malaysia.
Mean monthly temperatures at Zomba range between 18 and 25°C (Brakefield & Reitsma, 1991;Windig et al., 1994; Figure 2) and between 24 and 29°C at Sanguem (Figure 2). Both these locations are highly seasonal in terms of precipitation (Figure 2). At the third, the aseasonal Genting, monthly mean temperatures only vary between ~25 and 26°C and precipitation is received throughout the year ( Figure 2). Although environmental predictability is the core prerequisite for maintaining adaptive plasticity, few studies formally quantify the degree of predictability in relation to the evolution of phenotypic plasticity (Burgess & Marshall, 2014). Here, we have quantified the degree of seasonality and predictability of precipitation and temperature, and their synchronicity, by using wavelet analysis and Colwell's indices.
Based on the differences in environmental variation across the three study localities, we posit two hypotheses. Firstly, in Zomba and Sanguem, irrespective of variation in the temperature, strong seasonal precipitation will restrict host plant availability to the wet season. Thus, species are expected to synchronize their reproduction according to the seasons, with females breeding in the wet season, but arresting reproduction in the dry season (see Braby, 1995;Halali et al., 2020). In contrast, species from Genting will reproduce throughout the year as host plants are continuously available due to frequent rains. Secondly, the range of annual temperatures in Sanguem and Genting only represents the upper extreme of the temperature used in the common garden experiments to examine reaction norms and which always result in WSFs (e.g. van Bergen et al., 2017;Brakefield & Reitsma, 1991). Thus, the occurrence of distinct seasonal forms at Sanguem and Genting will indicate either that species have evolved a new polyphenic threshold for temperature or that they use cues other than temperature.

| Quantifying environmental predictability
We carried out wavelet analysis using the R package WaveletComp Ver 1.1 (Rösch & Schmidbauer, 2018a) to quantify the predictability of annual changes in precipitation, temperature and to investigate their synchronicity and phase differences for all three locations.
Wavelet analysis estimates frequency spectra at each time point in a time series revealing any abrupt changes in frequency through time (Burgess & Marshall, 2014;Cazelles et al., 2008Cazelles et al., , 2014Rösch & Schmidbauer, 2018b;Tonkin et al., 2017). Thus, a wavelet-based approach is superior to other widely used methods such as Fourier transform, which is well-suited for periodic time series data but it cannot detect abrupt changes in frequency over time (Cazelles et al., 2008). Furthermore, compared to Fourier transform, wavelet analysis is scale-independent and thus allows direct comparison of the predictability of multiple variables in a time series.
Firstly, we used wavelet analysis to characterize whether precipitation and temperature values recur during a defined period, for example, every 6-or 12-month period, for all three locations. Regularity in the events over a time series would suggest that the event is predictable. We extracted long-term data for monthly mean F I G U R E 2 Mean monthly precipitation (blue bars) and temperature (red line) for the three sampling locations. Climate data were extracted for the period from January 1901 to December 2019 as available from the Climatic Research Unit database precipitation and temperature of all three locations from January 1901 to December 2019 available from the Climatic Research Unit database (CRU) at ~55 × 55 sq. km. resolution (CRU TS Ver. 4.04; Harris et al., 2020). We visualized periodicity for both variables using the wavelet power spectrum (function wt.image). Regions of the significant period (p value < 0.05) were determined based on 100 simulations. To allow the comparison of both variables across all three sites on a similar scale, we first explored which combination of variable and location had the maximum power and then used the value, which is slightly higher than the square root of maximum power and 0.5 as an exponent (see Rösch & Schmidbauer, 2018b).
This approach allows making the colour gradient in spectral plots more explicit (Rösch & Schmidbauer, 2018b). We further plotted average wavelet power (function: wt.avg) for all three locations to determine the strength of periodicity in precipitation and temperature to a particular time period in a time series. Secondly, we carried out cross-wavelet analysis (function: analyse.coherency), which allows drawing conclusions about the nature of the synchronicity of two variables in a time series (e.g. precipitation and temperature), that is, whether they are in or out of phase. When variables are out of phase, this method also enables detailed quantification of the timing and nature of the phase difference by identifying which variable is leading and which is lagging with regards to the other (see Cazelles et al. (2008) and Rösch and Schmidbauer (2018b) for theoretical details and example R codes). We visualized coherency using wavelet power spectrum with additional information such as arrows on the plot showing details about the patterns with regards to phase difference. For example, arrows pointing straight to the right indicate that both series are perfectly in phase and arrows pointing straight left indicates that the series as completely out of phase. When the arrows move (rotate) between these two extremes, they will include a component pointing either upwards or downwards. This component indicates whether the first series (i.e. temperature) is leading (upwards) or lagging (downwards) on the second series (i.e. precipitation). We calculated the phase difference between precipitation and temperature using example R codes from Rösch and Schmidbauer (2018b).
Constancy measures the extent to which the environment remains similar for all months across years, whereas contingency measures the repeatability of seasonal patterns or the extent to which differences between months are similar across years (Colwell, 1974;Tonkin et al., 2017). Thus, contingency is a good measure of seasonality. Predictability is simply the sum of constancy and contingency; thus, higher predictability can be either due to high constancy and low contingency or vice versa. We, therefore, calculated the ratio of contingency and predictability to determine the extent to which seasonality contributes to predictability. We calculated all three indices using a time series data in the R package hydrostats Ver 0.2.7 (Bond, 2019). We used the default "transform" method with 11 breakpoints for binning data which is a required step for calculating Colwell's indices (Colwell, 1974).

| Long-term field sampling
Mycalesina butterflies were collected using traps baited with fermented banana. In Genting, we used collected samples (~8 samplings per month) from January 2012 to January 2013 in the previous study (van Bergen et al., 2016)  mnasicles and My. orseis are associated with forested habitats (Eliot et al., 1992). In Sanguem, three species (Mycalesis mineus, Mycalesis visala, M. leda) were abundant. Both My. visala and My. mineus occupy a range of habitats such as forests, scrubs, urban gardens and similar habitats (Kunte, 2000 Windig et al., 1994). Me. leda is the only species which occurred at both seasonal locations and is a widely distributed generalist species occupying a range of habitats from rainforests to human-dominated landscapes such as urban gardens (Kunte, 2000;Larsen, 2005).  et al., 2020). However, no similar trend was observed in Sanguem, probably due to lower sample size during the transition period, and hence all samples from the dry season were used in the analysis. We used samples from all wet season months in Zomba and Sanguem to calculate the proportion of reproductively active females in this season. Finally, as there is no distinct dry season in Genting, we visually (scatter plots) investigated whether there was any pattern of females undergoing diapause during any part of the year using samples from both sampling periods (see Figure 4).

| Measuring phenotypic traits
We carefully separated the wings from the thorax with micro-scissors and photographed using a Leica DFC495 camera coupled with a Leica M125 stereomicroscope (for Zomba and Genting samples) with an inbuilt scale or using a Nikon D5600 with Tamron 90 mm macro lens (for Sanguem samples) and graph paper as a scale. From these images, we measured the area of the CuA1 eyespot (using the yellow ring as the outer element of the eyespot) from the left or right least damaged ventral hindwing, and a proxy hindwing area using three landmarks defined by wing veins and then calculating the area of the enclosed triangle (see van Bergen et al., 2017, Figure S1). Relative eyespot size area was then calculated as the ratio of absolute eyespot area and proxy hindwing area. All images were analysed in ImageJ software (Schneider et al., 2012) using a custom-written macro. Re-measuring a subset of 50 specimens on different dates indicated high repeatability (R 2 = 0.99) for both eyespot and wing area.

| Statistical analyses
Apart from Me. leda, no species were sampled in more than one community. We, therefore, restricted our statistical analyses to within communities. We did not correct for phylogeny due to low sample size (total species = 10) and because phylogenetic comparative analyses focus on species means, whereas we are interested mainly in variation within species.
Many species of Mycalesina butterflies show a characteristic bimodal distribution in the eyespot size, with large spots in the WSF and small spots in the DSF (Brakefield & Larsen, 1984;Windig et al., 1994; Figure 1). We first explored the distribution of relative eyespot area (eyespot area/hind wing area) by plotting histograms for each species. We then classified our samples into WSF or DSF by setting a threshold value when eyespot area showed a strong bimodal distribution to separate the two peaks, with WSF and DSF above and below, respectively ( Figure 5). However, this procedure proved possible only for species from Sanguem and Zomba and not for those from Genting where distinct WSF and DSF forms were absent ( Figure 5, see below). Also, specimens from seasonal locations which were missing eyespots due to wing damage, we classified them into seasonal forms based on overall differences in the wing pattern (e.g. van Bergen et al., 2017).
We further explored seasonal variation in eyespot size using generalized additive models (GAMs) with Gaussian distribution separately for each location. GAM allows fitting a smoothing function to detect nonlinearity in the response variable (Pederson et al., 2019). The degree of nonlinearity is described by the effective degrees of freedom (EDF) with higher values of EDF indicating a greater nonlinearity, whereas the value of 1 suggests a linear trend (e.g. Hunsicker et al., 2016).
We included date (by first converting to number using a function in Microsoft Excel and then standardized to the mean of zero and unit standard deviation) and species as predictors, and relative eyespot size as a response variable. We fitted a smoothing term to the date and included species as an interaction term. We did not include sex as a predictor because the sample size for females was generally low and was entirely missing for some months in Sanguem and Genting. For Zomba, where we had a comparatively good representation of both sexes, a model including sex as a predictor had a marginally significant effect (p = 0.0271) and EDF values were similar to males (Table S1). GAM models were fitted using mgcv R package Ver 1.8.31 (Wood, 2019). All statistical analyses were performed in R Ver. 3.6.1 (R Core Team, 2019).

| Predictability of temperature and precipitation, and their synchronicity
Wavelet power spectra and average power plots for precipitation and temperature in Zomba showed strong periodicity indicating the presence of a single dominant frequency recurring every 12 months (Figures 2 and 3; Figure S2 and S3). In Sanguem, similar to Zomba, precipitation showed a strong periodicity over a 12-month period, whilst for temperature, two significant regions at 6-and 12month periods suggest that there are two prominent peaks which recur in every 12 months throughout the time series (Figures 2 and   3; Figure S2, S3). In Genting, both precipitation and temperature showed lower periodicity and had a lower average wavelet power compared to seasonal locations (Figure 3). Similarly, Colwell's indices suggested that precipitation and temperature in Zomba and Sanguem showed lower constancy and higher contingency or seasonality, but vice versa in Genting (Table 1).
Wavelet power spectra and average power plots for crosswavelet analysis indicated that both climatic variables were in synchrony and showed higher average power at the 12-month period ( Figure 3). Calculating the phase difference further showed that a rise in temperature always precedes the onset of rainfall, on average by 1.02 (range: 0.6-1.36) months in Zomba and 2.35 (range: 1.57-3.21) months in Sanguem. In addition, in Sanguem, the second peak (6-month period) was almost completely out of phase with precipitation ( Figure 3). Moreover, there was a trend of an increasing phase difference between temperature and precipitation at a 12-month period in Sanguem, suggesting that temperature has been rising substantially earlier in relation to the onset of rainfall in recent years; such an increase is not shown in Zomba ( Figure S5). In Genting, both variables were highly out of phase, with precipitation typically leading temperature by more than 3 months and they also showed a lower average power (Figure 3).

| Reproductive strategies in the dry season
In Zomba, combining data from three annual cycles, three previously studied species each exhibited a different reproductive strategy (see Halali et al., 2020). Briefly, B. campina reproduced throughout the year with 61% females in the mid-dry season having spermatophores and yolk-filled eggs, whereas less than 22% females at this time had spermatophores and eggs in B. safitza. Interestingly, in B. ena, 88% of females had spermatophores, but only 28% had eggs. An additional species from Zomba (Me. leda), which was not surveyed in the previous study and three species from Sanguem showed a similar trend with only <17% females having spermatophores and eggs in the dry season ( Figure 4). The only exception was My. mineus, where 41% of F I G U R E 3 Wavelet power spectra (three left panels) and average wavelet power plots (right panel) depicting periodicity of precipitation, temperature and synchronicity of both climate variables at each location. In wavelet power spectra plots, the strength of periodicity is indicated by the extent to which warmer colours are distributed over a particular time period. The black solid lines and white contour lines in spectral plots show regions of power significant at the 5% level based on 100 simulations. In average wavelet power plots, the strength of periodicity over a particular period is indicated by a higher power and the colour of each curve corresponds to the location. In the bottom panel, arrows in spectral plots indicate the phase difference between precipitation and temperature (see Methods) females had spermatophores but only 14% having eggs (Figure 4).
Moreover, in all species (except one) from Sanguem and Zomba,  (Figure 4). In all species, the presence of mature eggs was correlated with that of spermatophores in >90% cases.

| Seasonal changes in eyespot size
All species in Zomba and Sanguem showed a strong bimodal distribution for eyespot size with eyespots being large in the wet season and small or absent in the dry season ( Figure 5). Furthermore, EDF values from GAM analysis indicated strong nonlinearity in the eyespot size across sampling period (Table 2). In Genting, eyespot size exhibited unimodal distribution ( Figure 5), although EDF values suggested that the eyespot size of one of the species, My. intermedia, exhibited slight nonlinearity ( Figure S6).

| D ISCUSS I ON
Seasonal rainfall, which causes a strong temporal skew in the resource distribution, is one of the major drivers of the evolution of life-history strategies in animals (Halali et al., 2020;Tauber et al., 1986;Tonkin et al., 2017;Varpe, 2017). To cope with alternating selection pressures favouring reproduction or survival in the wet and dry seasons, respectively, animals time their life-history events or express season-specific morphological and life-history traits (e.g. Brakefield & Reitsma, 1991;McQueen et al., 2019;Tauber et al., 1986). However, to achieve a tight environment-phenotype matching, annual precipitation should exhibit high predictability and some environmental cue(s) should reliably signal the forthcoming season (Tonkin et al., 2017;Varpe, 2017).
Irrespective of local variation in mean monthly temperature, seasonal precipitation in both Zomba and Sanguem is expected to cause similar changes in the environment. Annual precipitation showed high predictability in both Zomba and Sanguem. It was synchronized with temperature, which was leading precipitation by at least a month. In contrast, such a seasonal trend was absent in Genting and here both variables were also highly out of phase. Such a characteristic pattern of temperature rising before the rainfall and dropping towards the end of the wet season, at least in East Africa (e.g. Zomba), is suggested to signal the forthcoming wet-and dry season, respectively, and hence acts as a reliable cue (Brakefield & Reitsma, 1991;Kooi & Brakefield, 1999;Windig et al., 1994). Whether temperature acts as a reliable cue in other parts of the tropics remains to be tested (see below). Thus, high predictability of precipitation and temperature acting as a cue (at least in Zomba, see below) in seasonal locations likely maintains seasonal polyphenism in the wild. Although, it should be noted that the sampling period differed between locations, from 1995 to 1998 in Zomba but more recently in Sanguem and Genting.
However, it is unlikely that this affects our results. Moreover, instead of using coarse climate data, microhabitat level data, especially for temperature, would be an invaluable addition in capturing more subtle consequences due to climate change. It is also interesting to note that we found an increasing phase difference between the lead time of temperature and on the rising precipitation in Sanguem in recent years ( Figure S6). Such a trend raises important questions about temperature being a reliable cue in signalling the forthcoming wet season given that the lag between the cue and future selective environment is becoming longer and exceeding the typical generation time of the organisms. Q uantifying the reproductive status of females based on the presence of spermatophores and mature eggs in both seasonal locations suggested that, except for a single species, females synchronized their reproduction according to the seasons. Thus, females actively reproduced in the wet season but arrested their reproduction in the dry season. The species breaking the pattern is B. campina, which belongs to a group of species, which are commonly found in drier forests but not in completely open habitats (Larsen, 2005;Windig et al., 1994). A nondiapausing strategy in this species, and also those in Genting, is likely due to more continuous availability of host plants in forests (e.g. Braby, 1995;Halali et al., 2020).
All species from Zomba and Sanguem exhibited a strongly bimodal distribution for eyespot size, a characteristic of seasonal polyphenism in Mycalesina butterflies, with eyespots being large in the wet system, and small or absent in the dry season (Brakefield TA B L E 1 Colwell's indices for precipitation and temperature for all three locations. Higher constancy indicates higher environmental stability; higher contingency indicates higher repeatability of seasonal patterns; predictability is the sum of constancy and contingency and contingency/predictability indicates the extent to which seasonality contributes to predictability  , 1991). Eyespots can function as anti-predatory devices which deflect predator attacks in the wet season when butterflies are active amongst green herbage, whilst an inconspicuous wing pattern in the dry season aids in crypsis of inactive butterflies amongst dry leaf litter (Lyytinen et al., 2004;Prudic et al., 2015;Stevens, 2005). Thus, strong selection for environment-phenotype matching in a seasonal environment likely maintains polyphenism in eyespot size (Brakefield & Reitsma, 1991;Lyytinen et al., 2004). In Genting, the distribution for eyespot size was more unimodal suggesting only minor variation across the years. However, investigating whether species from Genting have indeed reduced or lost plasticity in eyespot size will be more robustly answered by carrying out common garden experiments in the future. Interestingly, one species from Genting (My. intermedia) tended to show slightly smaller eyespots during January-March.
Overall, our data support the theoretical prediction that adaptive plasticity will be favoured in a more seasonal environment, but will be reduced or lost in a more stable environment, perhaps because there is a cost for maintaining plasticity (DeWitt et al., 1998;Reed et al., 2010;Tufto, 2015). Whether such costs exist in Mycalesina butterflies is not clear, but theory suggests that maintaining the sensory system required to process environmental information can incur some costs (Auld et al., 2010;DeWitt et al., 1998). Furthermore, large eyespots may function more effectively on the green background (Lyytinen et al., 2004;Prudic et al., 2015), and since it rains continuously in Genting with the habitat remaining green throughout the year, species are expected to maintain larger eyespot size.
However, it is interesting to note that two species (T. janardana and C. mnascicles) at Genting showed relative eyespot size similar to those of the DSF in species from the seasonal locations ( Figure 5). It is possible that these species may rely on alternative anti-predatory strategies.
Although seasonality is likely to be the major driver of the evolution of polyphenism in these butterflies (Brakefield & Larsen, 1984 Brakefield & Reitsma, 1991;de Jong et al., 2010;Nokelainen et al., 2018;Roskam & Brakefield, 1996). Such strong seasonal variation in temperature and the presence of distinct seasons likely explains the presence of seasonal polyphenism in species from F I G U R E 4 Mating strategies of females of species at the three locations (trends for Bicyclus species from Zomba are described in the Results section and see Halali et al. (2020)

F I G U R E 5
Histograms (bins = 40) showing distributions of relative eyespot size for all species including both sexes across the three locations (orange = Zomba; blue = Sanguem; green = Genting). The black vertical line on each plot indicates the threshold used for classifying seasonal forms such that individuals with an eyespot size below and above the threshold were classified as dry and wet season forms, respectively. Species from Genting were not classified into seasonal forms as distinct bimodal peaks were absent A recent study on a population of M. mineus from Thiruvananthapuram, India, which has a similar climate to Sanguem, demonstrated that humidity also affects pupal colour, a trait known to be phenotypically plastic (Mayekar & Kodandaramaiah, 2017).
When larvae were reared at 27°C with either high (85%) or low (60%) relative humidity resembling humidity during wet and dry seasons, respectively, the former condition resulted in more green pupae, whereas the latter resulted in more brown pupae. Plasticity in pupal colouration is likely to be adaptive as green pupae provide better background matching in the wet season and brown colour in the dry season (Mayekar & Kodandaramaiah, 2017). Another study revealed that seasonal fluctuation in relative humidity rather than the temperature was a better predictor of the temporal distribution of seasonal forms in the polyphenic butterfly Luthrodes pandava from India (Tiple et al., 2009). These studies indicate the possibility that humidity might be a better cue in regions where the temperature fluctuates within a narrow range. Furthermore, rather than a single cue, but possibly an interaction between multiple cues may provide a better prediction of the seasonal change (e.g. Goehring & Oberhauser, 2002;Morehouse et al., 2013;Singh et al., 2020).
For example, high temperature and low humidity can result in deteriorating plant quality, and larvae of Mycalesina butterflies growing on low host plant quality tend to have smaller eyespots (Kooi et al., 1996;Singh et al., 2020). In nature, host plant quality generally declines towards the end of the wet season, which may result in the expression of DSF (Brakefield & Reitsma, 1991 Our study demonstrates the advantage of quantifying the extent of predictability, as well as the nature of correlations between potential cues, in temporal variation in climatic variables when making comparative analyses of phenotypic plasticity across different communities and locations. We conclude that for Mycalesina butterflies, seasonal precipitation, which results in distinct selective environments, and the predictability of precipitation are the primary factors influencing the evolution of polyphenism (e.g. Roskam & Brakefield, 1999).
Populations or species may secondarily evolve local adaptations for cue(s) use that depend on the local variation in the environment (e.g. Lindestad et al., 2019;Roskam & Brakefield, 1999). In addition, strong temperature-dependent expression of seasonal forms in Mycalesina butterflies as revealed by laboratory experiments may provide biased results since the temperature is often the only variable tested (but see Mayekar & Kodandaramaiah, 2017;Singh et al., 2020). Future reaction norm studies using multiple factors such as temperature, humidity or host plant quality in a full factorial design will be important in helping to understand the regulation of seasonal polyphenism in the wild. Cue use has special relevance in estimating the consequences of climate change as it is likely to affect their reliability in predicting the forthcoming environment under changing patterns of seasonality (Bonamour et al., 2019). Thus, understanding local variation in cue use will help in predicting future responses of communities of Mycalesina butterflies and other tropical polyphenic taxa to climate change.
F I G U R E 6 Phylogenetic relationships of Mycalesina species represented in the study. Note that currently there is no phylogeny that includes Me. leda together with all genera of Mycalesina butterflies but the species is likely to be ancestral to Mycalesina (see Chazot et al., 2019). The location of species is indicated by species colour codes and corresponding circles indicate whether species are polyphenic or non-polyphenic based on the presence/ absence of strong bimodal distribution for the eyespot area (see Figure 5)

CO N FLI C T S O F I NTE R E S T
Authors have no conflicts 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.13895.