Population biology and phenology of the colour polymorphic damselfly Ischnura elegans at its southern range limit in Cyprus

1. Geographically widespread species provide excellent opportunities to investigate how phenotypes change across large‐scale environmental gradients. Temperature is a fundamental environmental variable and an important determinant of insect fitness. However, field research is often geographically restricted, and typically concentrated in northern latitudes. Basic population biology and phenotypic clines in relation to temperature therefore remain poorly known across the entire geographic range, even in otherwise well‐studied taxa.


Introduction
The genetic structure and reproductive phenology of animal populations vary geographically, particularly in widely distributed species that occupy diverse habitats and encounter contrasting ecological conditions (Slatkin, 1987;Conover, 1992;Lenormand, 2002;Morrison & Hero, 2003). Quantifying such Correspondence: Erik Svensson, Evolutionary Ecology Unit, Department of Biology, Lund University, SE-223 62 Lund, Sweden. E-mail: erik.svensson@biol.lu.se geographic variation is necessary to understand how geographic differences in selection can shape genetic and phenotypic diversity. For example, latitude is one of the most intensely investigated geographic axes of genetic differentiation between animal populations (Palo et al., 2003;Fabian et al., 2012). Latitude is also associated with clinal variation in numerous phenotypic traits (Azevedo et al., 1998;Schmidt & Paaby, 2008;Parsons & Joern, 2014), the timing of breeding periods (Schmidt & Paaby, 2008;Mazaris et al., 2013;Parsons & Joern, 2014) as well as allele frequencies underlying discrete phenotypic morphs (Forsman & Shine, 1995;Phifer-Rixey et al., 2008;Gosden et al., 2011;Takahashi et al., 2011). A better understanding of how environmental factors such as temperature, seasonality, and humidity underpin observed latitudinal gradients can be used to capitalize on geography to infer ecological and evolutionary responses to ongoing climate change (De Frenne et al., 2013). Yet, to achieve this goal, an essential first step is to build up a comprehensive record of population features across diverse species and their entire distribution ranges.
Despite the ubiquity of geographic variation in mean phenotypes, allele frequencies, and phenologies, most research on many wide-ranging species have been limited to the areas inhabited by the majority of researchers, leading to a geographically restricted understanding of their ecology (Zuk, 2016). In terms of numbers of researchers and more favourable economic conditions for basic research, there has been an ongoing bias to carry out research in northwestern Europe and North America, whereas other regions of the world remain largely understudied. For example, the Common Bluetail, Ischnura elegans, is a pond damselfly species (family Coenagrionidae) with a broad distribution across the Palearctic. It occurs from Spain in the West, to Japan in the East, and from Sweden in the North to Iran in the South (Boudot & Salamun, 2015). Classic studies on the life-history and demography of I. elegans (Parr, 1965(Parr, , 1973Parr & Palmer, 1971;Parr & Parr, 1972;Van Noordwijk, 1978;Hinnekint, 1987;Anholt et al., 2001) have been all conducted in Western Europe and Scandinavia, thereby capturing only a small fraction of the environmental conditions this species encounters across its entire range. Similarly, the flight season, during which individuals may be reproductively active and subject to sexual selection, has been investigated in depth only in this same region (Parr, 1965(Parr, , 1970. Although I. elegans has been subject to many recent ecological and evolutionary studies under natural or seminatural conditions (Takahashi et al., 2014;Le Rouzic et al., 2015;Lancaster et al., 2017;Willink & Svensson, 2017;Willink et al., 2019;Svensson et al., 2020), there is a strong geographic bias in these previous studies, which limits our understanding of its reproductive and population biology, and how they vary across a larger geographic area.
During the last two decades, I. elegans has been subject to numerous studies focusing on the behavioural mechanisms behind, and the population consequences of, the maintenance of genetic variation (Van Gossum et al., 1999Svensson et al., 2005;Le Rouzic et al., 2015). This interest is motivated by the existence of a heritable and discrete female-limited colour polymorphism, which can be used as a visual genetic marker to study I. elegans (Sánchez-Guillén et al., 2005;Svensson et al., 2009). Previous studies have shown that the female colour polymorphism in I. elegans is maintained locally by negative frequency-dependent selection via male mating and pre-mating harassment Le Rouzic et al., 2015). Males disproportionately target their mating attempts towards a particular female morph as its frequency increases (Van Gossum et al., 2001). Excessive mating in turn reduces female fecundity (Gosden & Svensson, 2007), leading to a fitness advantage to rare females . In this and other damselfly species with similar female-colour polymorphisms, one of the female morphs has a colour pattern nearly identical to that of males (Robertson, 1985;Van Gossum et al., 2011). Such male-like females benefit from their visual resemblance to males, which acts to further reduce superfluous male-mating attempts, and they are therefore considered to be male mimics (Robertson, 1985;Hammers & Van Gossum, 2008;Iserbyt et al., 2011;Van Gossum et al., 2011). However, due to plasticity in male search image formation, the extent of such an advantage of male mimicry is also frequency-and probably density dependent (Hammers & Van Gossum, 2008;Iserbyt et al., 2011). These frequency-dependent male-female mating interactions will prevent the fixation of any single female morph locally. In fact, of over one hundred populations surveyed in Europe, none is monomorphic and only a few (∼10%) appear to be dimorphic, with one of the three female morphs being absent or extremely rare (Gosden et al., 2011). This suggests that such frequency-dependent selection is strong and operating in most populations. Moreover, in a decade-long longitudinal study of more than a dozen populations of I. elegans in southern Sweden, none lost any of the three morphs and all remained trimorphic during more than ten generations (Le Rouzic et al., 2015).
Despite the frequency-dependent processes that maintain the polymorphism locally, populations can vary quite markedly in morph composition across large geographic scales (Gosden et al., 2011). Importantly, female morphs differ in several phenotypic traits, besides colouration, that are ecologically relevant, namely resistance and tolerance against parasitic mites (Willink & Svensson, 2017), cold tolerance (Lancaster et al., 2017), larval developmental time  and temperature sensitivity of adult sexual maturation and colour development . Therefore, ecological selective agents, other than male mating harassment, such as climatic factors, parasite prevalence, and virulence, are also likely to shape geographic variation in morph frequencies across the large range of I. elegans. For instance, the equilibrium frequency of the more cold-tolerant morph could be higher in more northerly populations, with cooler temperatures during the flight season . Such ecological agents shaping regional variation in morph frequencies are also likely to operate in tandem with sexual conflict dynamics. Our knowledge about how these different selection processes interact, reinforce or oppose each other can be increased with broader geographic studies of the population biology of I. elegans, to include the full range of the ecological conditions that this widely distributed species encounters across its range. In particular, in-depth studies at the southern range limit of this widespread species should complement recent studies on phenotypic and genomic variation at its northern range limit in Sweden (Lancaster et al., 2016(Lancaster et al., , 2017Dudaniec et al., 2018). In addition, the environmental and genetic factors that determine the location of species range limits is a question that is of central interest to all ecologists, evolutionary biologists and conservation biologists (Kirkpatrick & Barton, 1997;Bridle & Vines, 2007).
Here, we report the results from surveys of breeding populations of I. elegans throughout Cyprus, at the southern range limit of this species in Europe. We recorded female morph frequencies using standardised field sampling procedures based on previous studies Le Rouzic et al., 2015). These data add to a growing database of standardised population samplings of I. elegans in Europe where morph frequencies have been quantified in a similar fashion . Moreover, Cyprus has a Mediterranean climate, that differs strongly from the climate in Sweden and northwestern Europe (Cornes et al., 2018), where most previous studies on natural populations have been conducted.
We also report on the flight season of I. elegans on Cyprus and compare it to the flight season in other parts of Europe. Flight season data were derived from the Cyprus Dragonfly Study Group (CDSG) database, from records collected between January 2013 to December 2019. Population densities and the strength of mating interactions are likely to vary seasonally and geographically and are likely to be influenced by the overall length of the reproductive season. Consequently, information on reproductive phenology is important for understanding the interaction between natural and sexual selection in this and other animal species. Furthermore, by influencing the strength and opportunity for intersexual interactions, breeding phenology may influence the fitness costs and benefits of the alternative reproductive and ecological strategies in the different morphs (Machado et al., 2016), thereby shaping geographic variation in the equilibrium morph frequencies among populations.
Finally, we present data on body size in both sexes and all three female morphs from Cyprus. Adult body size is a fundamental and ecologically important phenotypic trait that is almost invariably a target of natural and sexual selection in dragonflies and damselflies (Waller & Svensson, 2017). Our study therefore provides the groundwork necessary to establish I. elegans as a model organism for ecological studies across wider biogeographic scales.

Study species
Colour development in females of I. elegans is controlled by a single locus or a set of tightly linked loci, with sex-limited expression and three alleles in a dominance hierarchy (Sánchez-Guillén et al., 2005;Svensson et al., 2009). Thus, in females, the six possible genotypes give rise to three discernible colour morphs, whereas males are always monomorphic. One of the female morphs, traditionally known as Androchrome and hereafter referred to as A-females, develops a colour pattern very similar to that of males, with green-blue antehumeral stripes and a bright blue patch on the eighth abdominal segment (Cordero et al., 1998;Svensson et al., 2009;Fig. 1a). The other two female morphs, denoted as Infuscans and Infuscans-obsoleta (and abbreviated as Iand O-females respectively), are strikingly different from males when sexually mature (Cordero et al., 1998;Svensson et al., 2009;Fig. 1b,c). I-females in their final, sexually mature colouration have olive-green to brown antehumeral stripes (Fig. 1b). O-females lack stripes but exhibit a uniform pink thoracic colour (Fig. 1c), that turns bronze-brown upon sexual maturity. Both I-and O-females express the bright blue abdominal patch of males and A-females early after their emergence from the larval state . However, these two female morphs conceal the blue colouration during adult development, so that unlike males and A-females, sexually mature I-and O-females exhibit a dark abdominal colouration dorsally (Cordero et al., 1998;Svensson et al., 2009;Willink et al., 2019). Immature A-and I-females have similar violet markings on the thorax, and their striking colour differences described above develop with the onset of sexual maturation (Willink et al., 2019. The Common Bluetail, I. elegans, breeds in diverse habitats throughout Europe. It is found in ponds, lake shores, and slow-flowing streams, primarily in open landscapes and up to an elevation of 1600 m.a.s.l. (Boudot & Salamun, 2015). I. elegans is considered a widespread, non-threatened species within its local range and is highly tolerant to human disturbance of aquatic habitats (Harabiš & Doln, 2012). Voltinism, the number of generations in a year, is inversely correlated with latitude, ranging from larvae which take more than a year to emerge onto land in Scotland, to multiple generations in the same season in Southern Europe (reviewed in Corbet et al., 2006). Latitudinal variation in voltinism is in turn contingent on variation in the length of the reproductive season. Reproductive phenology has not been systematically investigated in I. elegans, yet in other widely distributed species of dragonflies and damselflies breeding seasons are longer at lower latitudes (reviewed in Corbet, 1999).  Le Rouzic et al., 2015). In short, a single researcher walked a transect along suitable breeding habitat in each locality, collecting individuals with a hand net upon encounter. Collected individuals were temporarily housed in plastic mesh cages if single and in plastic cups, if collected while mating. After the sampling period (mean = 32.58 min; SD = 19.78), we recorded the sex and mating status of all collected individuals, and for females, the colour morph and developmental colour phase were also recorded. Females were classified to the immature colour phase if they displayed a bright blue colour patch and, in the case of A-and I-females, if they had violet markings in the thorax. Females with their final sexually mature colour pattern were classified to the mature colour phase.
Individuals were subsequently released at or near (within a few hundred meters) their site of collection. Temporarily caging individuals prevented us from double-counting any damselflies and allowed for more accurate estimates of female morph frequencies. A few sampled individuals (mean = 2.42; SD = 1.39; per sex and female morph) in each of the localities with a higher abundance of I. elegans (n = 26) were used for body size measurements prior to their release. We used a hand-held electronic calliper to measure the total length (excluding genital appendages), abdomen length, and hind wing length of these individuals. The CDSG is a citizen-science group founded in late 2012. Since 1st January 2013 its members have been monitoring dragonfly and damselfly communities along defined transects at over 50 sites, selected to give broad geographic coverage, and include all habitat types and known species in the island. For each species at each site, the number of single individuals is counted, where possible recording males and females separately. The numbers of pairs in tandem, pairs in copula, and females ovipositing are also recorded. These records are entered into the CDSG database along with other incidental records and those from visitors and non-CDSG members. The CDSG aims to monitor each site monthly or twice-monthly, yet as this is a citizen-science effort, data completeness varies among sites. Although several of the initially selected sites were lost for a variety of reasons, such as drying up in years of low rainfall or becoming overgrown, particularly with the Mediterranean Reed (Arundo donax), new sites were found and included in the monitoring. Between 1 January and 31 December 2019, a total of 7877 site visits had been made. Here, we use these records to estimate the country-wide flight period, following the standardised procedure in Boudot and Kalkman (2015) for European Odonata. To investigate seasonal patterns in the flight period of I. elegans across localities, we used the sites where monthly data is available for at least 36 months (up until 25 June 2019). Both phenology and morph frequency surveys were conducted during the activity hours of I. elegans between 08.00 and 15.00 avoiding periods of rain and strong winds.

Statistical analysis
We used generalised mixed-effect models fitted by MCMC to estimate the frequencies and mating rates of the three female morphs of I. elegans in Cyprus. The analyses were conducted using the package MCMCglmm (Hadfield, 2010) in R (R Core Team, 2018). To estimate island-wide morph frequencies, we used a multinomial model accounting for the random effects of the sampling locality during each mating season. The residual variance (not identified in a multinomial model) was fixed to 1 and the residual covariance was set to 1/3 * (I + J), where I and J are two-dimensional identity and unit matrices (Hadfield, 2018). This means that we start with a priori belief that every sampled female has a similar probability of being of any of the three morphs. We used parameter expanded priors on the morph-specific between-locality variances to improve mixing of the variance components (Hadfield, 2018). These priors come from a non-central scaled F-distribution, and we used a scale of 25 following Gelman (2006).
To estimate morph-specific mating rates we used another categorical model with mating as a binary (0 or 1) response variable and morph as the fixed predictor. Here, we also fixed the residual variance to 1, but there was no prior on the residual covariance because the binary response can be reduced to a single dimension. Unlike the previous model, we allowed for different and independent variance components for the three sampling seasons (2017-2019) in each locality. This way, we accommodated variation in mating rates due to factors such as weather conditions, which could vary between sampling seasons and in different ways in different localities. Also, we did not specify a random categorical interaction and instead assumed that population and season effects on mating rates were similar across female morphs. Finally, we used a normally distributed fixed effect prior with mean equal to 0 and standard deviation equal to 1 + π 2/3 , as this prior is relatively flat on the probability scale when using a logit link (Hadfield, 2018).
We present the results of analyses including only those localities in which at least 10 females were sampled. For estimation of morph frequencies, we included females of both developmental colour phases, whereas for the estimation of mating rates we considered only females in their final, sexually mature colour phase. Both models were run for 2 250 000 iterations thinning every 1000th and with a burn-in of 250 000. We report posterior means (PM) and 95% highest posterior density interval (HPD) of parameter estimates. As significance tests of mating rate differences between morphs, we report the proportion of the posterior estimates in which the female morph with the lower mean mating rate has a higher mating rate than the other morph. We consider a difference in mating rate significant if this proportion was below 5%. We performed MCMC diagnostics for both morph-frequency and mating rate analyses using the R package coda (Plummer et al., 2006). We visually evaluated the stationarity of posterior distributions and convergence of independent runs. We also evaluated convergence using the Gelman-Rubin diagnostic and report effective sample sizes and autocorrelation between posterior samples for fixed effects (Figures S1-S2).
We estimated sex-and morph-specific body size (total length, abdomen length and hind wing length) and compared it to the range reported for Europe and the U.K. (Dijkstra & Lewington, 2006). We used a Gaussian model for each size measurement, also including random effects of the different localities and seasons sampled. In these models, we used a normal prior with mean of 0 and large variance (10e+8) on the size estimates and inverse Wishart priors with a low degree of belief on the residual and random variance components. The MCMC sampling strategy for these models was the same as for the morph frequency and mating rate analyses. MCMC performance was also diagnosed as mentioned above ( Figure S3). As significance tests, we report the proportion of the posterior distribution of body size estimates that fall within the European/U.K. range, for males and for each of the female morphs. We considered a difference in size significant if this proportion was below 5%.
We determined the flight period of I. elegans in Cyprus from the number of records (2793 in total) entered in the CDSG database for the period from January 2013 to December 2019. We followed the protocol used in Boudot and Kalkman (2015), except we grouped records half monthly (i.e. 1st to 15th and 16th to 28th-31st of each month) rather than in 10-day periods to reflect the monitoring schedule. Since the earliest and latest sightings often refer to unusual events, we define the start and end of the flight season as the first half month in which one or 99 percent of the records respectively have been made (Boudot & Kalkman, 2015). Here, a record is defined as an observation of I. elegans during sampling by the CDSG, independently of the number of individuals observed. Boudot and Kalkman (2015) normally define the main flight season as the first and last half month period in which 10 percent or more of the total records occur. However, I. elegans has a long flight season in Cyprus and the 10 percent level is not achieved in any half month. Thus, we lowered the threshold for the main flight season to five percent of the total records, as Boudot and Kalkman (2015) did with other species with a long flight season in Europe.
To investigate patterns of seasonality in the flight period of I. elegans in Cyprus, we focused on the sites with the most complete data over the period between January 2013 and June 2019, as well as the island-wide pattern. We selected localities with data for at least 36 months, and combined two pairs of such localities in very close proximity ('Morphou area 1' and 'Morphou area 2' into 'Morphou area' and 'Diarizos River 1' and 'Diarizos River 2' into Diarizos River 1-2, n = 8). However, as population surveys are conducted by volunteers, there were periods of missing data in all localities ( Figure S4). A total of 750 sampling events were conducted by the CDSG in these core localities during the study period. For each site, we obtained the average number of adults recorded per observer per sampling per month. We interpolated missing values in these time series using the na.approx function in the R package zoo (Zeileis & Grothendieck, 2005). We then obtained three measures of seasonality for each site and for the island-wide data using the package greenbrown (Forkel et al., 2013;Forkel et al., 2015). Here, the 'pgram' test for seasonality is positive if the maximum frequency of a periodiogram based on a fast fourier transformation equals 1 (i.e. one period every year). A positive 'acf' test indicates a minimum autocorrelation in a de-trended time series at a lag of 0.5 (i.e. the least correlation between surveys is observed at six-month intervals). Finally, the 'lm' test fits two linear models to the time series, one including just the trend and one including the trend and seasonal cycle. The test is positive if the second model is selected based on the Bayesian information criterion (BIC), that is, if inclusion of the seasonal cycle parameter substantially increases model fit. In all cases, these tests of seasonality yielded identical results if applied to either

Female-morph frequencies and mating rates
We visited a total of 45 localities, of which 19 had high enough female density to sample 10 or more females over the course of up to three mating seasons (Tables S1 and S2). In these populations, we sampled a total of 1682 individuals of which 619 were females. A-females were rare, if at all observed in these localities (Fig. 2a,b). The island-wide PM frequency of A-females was only 5.45% [95% HPD (0.038, 0.074)], whereas I-females [PM = 0.462; 95% HPD (0.398, 0.523)] and O-females [PM = 0.484; 95% HPD (0.419, 0.550)] occurred at similarly high frequencies (Fig. 2a,b). This pattern was mirrored across individual localities ( Figure S5). The two common female morphs also had higher mating probabilities than A-females (Table 1; Fig. 2c). The PM mating probability Table 1. Statistical contrasts of mating rates between female morphs of I. elegans in Cyprus. 'PM' represents the posterior mean of the difference in mating rate between a pair of morphs, '95% HPD interval' is the 95% highest posterior density interval of the difference estimate, and 'PMCMC' is proportion of the posterior distribution in which the female morph with the lower mean mating rate has a higher mating rate than the other morph.

Body size
We obtained body size data for 204 individuals of I. elegans in Cyprus (N males = 90; N A-females = 10; N I-females = 57; N O-females = 47). Both males and females in Cyprus were smaller in total length than the minimum reported for Europe and the U.K. (all PMCMC <0.02; Fig. 3a). I-and O-females were both larger than males (both PMCMC ≤0.001; Fig. 3a) and similar to each other (PMCMC = 0.192; Fig. 3a), whereas A-females were of intermediate size, overlapping in total length with males, I-females and O-females (all PMCMC >0.05; Fig. 3a). Despite this overall difference in size, the abdomen length of I. elegans in Cyprus was within the known European/U.K. range, although within the lower half of the range (Fig. 3b). Male abdomens were also shorter than those of I-and O-females (both PMCMC < 0.001; Fig. 3b), while A-females were intermediate (all PMCMC > 0.05; Fig. 3b). Females of all morphs had longer wings than males (all PMCMC < 0.001), and were within the reported European/U.K. range (Fig. 3c). In contrast, the posterior distribution of male hind wing lengths was centred around the range minimum for Europe and the U.K. (PMCMC = 0.538; Fig. 3c).

Flight season
During the period between January 2013 and December 2019 members of the CDSG had made a total of 7877 site visits to around 703 localities island-wide. I. elegans was recorded during 2788 of these visits with an adult count of 56169 individuals. I. elegans is by far the most common zygopteran on Cyprus and has been recorded from 216 localities island-wide (Fig. 1d), from sea level up to the Prodromos reservoir at 1600 m.a.s.l., which is the highest water-body monitored by the CDSG.
The flight season of I. elegans in Cyprus starts in the second half of February and continues to the second half of November (Fig. 4a). The main flight season is from the second half of March to the first half of September (Fig. 4a). This is substantially longer than that reported for the more northerly regions (Fig. 4a). There was evidence of seasonality in the island-wide occurrence data according to the periodicity ('pgram') and autocorrelation ('acf') criteria, but not according to the linear model ('lm') criterion (Table 2, Fig. 4b), which compares the fit of models with and without a seasonal cycle via BIC. Evidence for seasonality varied among localities, with some localities showing a seasonal pattern under periodicity or both periodicity and autocorrelation criteria (Table 2). However, the 'lm' criterion of seasonality was not met in any locality with long-term data (Table 2).

Discussion
I. elegans is a geographically widely distributed damselfly species that occurs across the Palearctic (Boudot & Salamun, 2015) and has become a well-known study system in evolutionary ecology. Population studies on I. elegans have addressed a variety of topics, including the age-, frequencyand density-dependent dynamics of male-mating harassment Van Gossum et al., 2011;Willink et al., 2019) and the maintenance of heritable morphs by balancing selection (Takahashi et al., 2014;Le Rouzic et al., 2015). Population studies with I. elegans have also been used to investigate phenotypic correlations between alternative mechanisms of defence against parasites (Willink & Svensson, 2017), and the role of environmental and social effects on temperature sensitivity and range expansion (Lancaster et al., 2017). These studies rely on female colour morphs as visual genetic markers of suites of correlated traits on which natural and sexual selection operate. However, the majority of such field studies have taken place near the northern limit of the distribution range of I. elegans. How these phenomena are influenced by large-scale environmental variation is therefore an important question, which could be addressed with comparative studies across the geographic range of I. elegans. As a first step on the road towards such a broader geographic scope of studies of I. elegans, here we present data Table 2. Seasonality in the flight activity of I. elegans across eight localities and for the island-wide data. Seasonality was tested according to three criteria (see Methods): maximum frequency of a periodiogram ('pgram'), the minimum autocorrelation of a de-trended time series ('acf'), and a model comparison with and without the seasonal component ('lm'). 1 = seasonality detected; 0 = seasonality not detected. on the basic population biology and phenology of this trimorphic species in Cyprus. This is the southernmost region where populations of this widespread species have been systematically studied to date. Three main features distinguish breeding populations in Cyprus from those in Northern Europe. Firstly, male-mimicking A-females are the minority morph in Cyprus, occurring only at ∼5% frequency, while I-and O-females occur at similarly high frequencies (Fig. 2a,b). Such a low frequency of A-females is striking, given that this may be the most visually conspicuous morph to human observers. A-females display a bright blue colouration throughout their adult life, which may increase their detectability compared to the other morphs, which develop a darker and duller colour pattern during sexual maturation (Henze et al., 2019;Willink et al., 2020). In contrast, the frequency of A-females increases with latitude, with females in Southern Sweden typically being composed of 60-80% male mimics (Gosden et al., 2011;Le Rouzic et al., 2015). The increasing frequency of A-females in northern Europe is likely due to their developmental advantages at cooler temperatures, whereby A-females enjoy higher pre-reproductive survival and faster sexual maturation and colour development . The developmental success of I-and O-females after emergence from the last nymphal stage is in contrast more sensitive to temperature, and these two morphs are therefore expected to benefit more from the warmer Mediterranean climate of Cyprus .
The low island-wide frequency of A-females might imply that some populations in Cyprus are effectively dimorphic, with only I-and O-females. In fact, a previous survey across continental Europe suggested that about 10% of populations were dimorphic, although they all included A-females (Gosden et al., 2011). The loss of a female morph may cause increased pre-mating and mating harassment in the other two morphs, as males would have fewer competing targets while forming a search image of potential mates (Dukas & Kamil, 2001). However, if male-mimicry is effective, the local absence of otherwise rare A-females should not dramatically alter male-mating harassment in the other two morphs, which would already account for most male-mating attempts under negative frequency-dependent selection. The local absence of A-females might also be temporary. In I. elegans, the A-allele is dominant over I-and O-alleles. While dominant alleles under negative frequency-dependent selection are more likely to be lost by genetic drift, they also re-invade populations more easily due to Haldane's Sieve, the expectation that a weakly advantageous mutation will increase more rapidly in frequency if dominant (Haldane, 1924;Pannell et al., 2005). Haldane's Sieve acting on migrants might thus contribute to the maintenance of the dominant A-allele at a regional scale (Roux & Pannell, 2019).
As in Cyprus, A-females are usually not the majority morph in the south of continental Europe, but typically I-females are more common than O-females (Gosden et al., 2011;Svensson et al., 2020). In the closely related species Ischnura genei, populations in the Mediterranean island of Sardinia also have relatively low frequencies of A-females, but O-females are generally more common than I-females (Sanmartín-Villar & Cordero-Rivera, 2016). Although several ecological differences between A-and both I-and O-females have previously been reported for I. elegans (Lancaster et al., 2017;Willink & Svensson, 2017;Svensson et al., 2020), ecological differences between I-and O-morphs have not been investigated to a similar extent, probably because O-females are so rare in northern Europe. The role of ecological mechanisms versus historical contingencies and 'island effects' that might shape morph-frequency variation within the Mediterranean region therefore remains an interesting question that should be addressed in the future, to get a better understanding of the ecological factors and evolutionary processes operating in this trimorphic system.
Secondly, adult individuals of I. elegans are relatively small in Cyprus (Fig. 3). This is particularly the case for males, which are also generally smaller than females , and when compared to the reported European range of total body length and hind wing length (Dijkstra & Lewington, 2006;Fig. 3b,c). This qualitative result is consistent with Bergman's rule, a pattern of increasing body size with decreasing temperature. Bergman's rule is supported in vertebrate endotherms (Ashton et al., 2000;Meiri & Dayan, 2003;Salewski & Watt, 2017; but see Riemer et al., 2018), but is not generally supported across diverse clades of insects and other ectotherms (Mousseau, 1997;Blanckenhorn & Demont, 2004;Adams & Church, 2008;Shelomi, 2012;Wonglersak et al., 2020). One mechanistic explanation for Bergman's rule that is applicable for I. elegans, and odonates in general, is that developmental rate (i.e. cell division and differentiation) increases more rapidly with temperature than does metabolism (Blanckenhorn & Demont, 2004). Therefore, higher developmental temperature at more southern latitudes should result in faster maturation at a smaller body size. Because growth in odonates occurs only during the aquatic larval stage, differences in developmental temperature should affect both the duration of the growing period (i.e. voltinism) and the adult body size (Johansson, 2003;De Block et al., 2008;Hassall et al., 2014). Our study did not address whether Bergman's rule is met throughout the distribution range of I. elegans. However, a decreasing number of generations per year at higher latitudes has been previously reported (Corbet et al., 2006), and the extended flight season of I. elegans in Cyprus (see below) also suggests that populations on this island are multivoltine.
A-females in Cyprus were more male-like in size than either I-or O-females. This is consistent with previous analyses in Swedish populations, showing that A-females are more male-like in shape than the other morphs ). Recent studies in female-polymorphic insects, including the widespread tropical and subtropical I. senegalensis (a congener of I. elegans that has been intensively studied in Japan), suggest that the development of male-coloured females is more masculinised, compared to the development of alternative female morphs (Takahashi et al., 2019). These developmental differences between morphs may be caused by alternative splicing and differential expression patterns of the regulatory gene doublesex (Takahashi et al., 2019(Takahashi et al., , 2020, which also underlies the development of somatic sex differences across many insect taxa (Kopp, 2012). To date, there is no direct evidence of a male-like expression pattern of dsx in A-females of I. elegans. However, the locus or set of tightly linked loci that govern colour morph development in I. elegans seems to have pleiotropic effects during colour development and differentiation of the female morphs . Such widespread pleiotropy may also impact the rate and duration of larval development, in turn generating size differences between morphs. In southern Sweden, the larval developmental period is shorter in the offspring of O-females, suggesting this allele is associated with a faster developmental rate . In contrast, the A-allele may be associated with a slower growth rate, as A-females have a similar developmental period as I-females but become mature at a smaller size . Whether these developmental differences occur at warmer temperatures has not been investigated.
Finally, the flight season, during which adult damselflies emerge and can potentially mate, is considerably longer in Cyprus than in any other location with comparable data (Fig. 4a; Boudot & Salamun, 2015). This is consistent with overall fast development and multivoltinism caused by warm to mild temperatures throughout the year. Flight activity in Cyprus spanned more than 9 months (Fig. 4a), resulting in partial evidence for a seasonal pattern in the island-wide dataset (Table 2; Fig. 4b). Some core localities with long-term data, had even less support for seasonality in flight activity (Table 2). Variation in seasonality among localities could be explained by regional differences in the length of growth periods and breeding seasons, for instance due to an altitudinal temperature gradient. Although a test of this hypothesis would entail estimating seasonal strength from a larger number of localities with long-term data, such a parallel between latitudinal and altitudinal gradients in the seasonality of (potential) mating activity is known for other ectotherms (Morrison & Hero, 2003). Nevertheless, the marked contrast between the flight season of I. elegans in Cyprus and more Northern European sites (Fig. 4a) suggests a pervasive role of temperature influencing life-history evolution in this widely distributed damselfly.
In conclusion, species with broad distribution ranges, such as I. elegans, provide excellent opportunities to investigate how large-scale climatic variation shapes the phenotypic outcomes of selection driven by local interspecific and intraspecific interactions. However, for many widespread species, local field studies tend to closely match the geographic distribution of the scientists who study them, which in turn reflects economic factors and research traditions, rather than strict biological considerations, and can lead to biases in the perception of which ecological factors are most important (Zuk, 2016). Here, we have studied the population biology of the widely distributed damselfly I. elegans in Cyprus, which is the southern range limit of this species in Europe. Populations of I. elegans in Cyprus are distinguished by the rarity of male-mimicking females, reduced body size, and a long flight season. These ecological differences from northern Europe, where the majority of field studies have been conducted, underscore the importance of broadening the geographic scope of field studies in I. elegans and many other widespread organisms.

Data availability statement
The data that support the findings of this study and all R code to reproduce the analyses are available online at https://doi.org/10 .6084/m9.figshare.13119980.v2.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. Traces for two independent runs of a model estimating female morph frequencies of the Common Bluetail damselfly (Ischnura elegans) in Cyprus. Traces are shown for two morphs as the third frequency can be obtained from the other two. Autocorrelation between posterior samples is <0.05, and effective sample sizes = 2000 in both cases. The multivariate potential scale reduction factor (mpsrf) of the Gelman-Rubin convergence diagnostic is 1.01. Figure S2. Traces for two independent runs of a model estimating female morph-specific mating rates in the Common Bluetail damselfly (Ischnura elegans) in Cyprus. Autocorrelation between posterior samples is <0.05, and effective sample sizes >1845 in all cases. The multivariate potential scale reduction factor (mpsrf) of the Gelman-Rubin convergence diagnostic is 1.01. Figure S3. Traces for two independent runs of models estimating body size in the Common Bluetail damselfly (Ischnura elegans) in Cyprus. Autocorrelation between posterior samples is <0.07, and effective sample sizes >1612 in all cases. The multivariate potential scale reduction factor (mpsrf) of the Gelman-Rubin convergence diagnostic is 1.01 in all three models.