Landscape context and plant population size affect morph frequencies in heterostylous Primula veris—Results of a nationwide citizen‐science campaign

Heterostyly is a genetically determined floral polymorphism of style length promoting outcrossing between individuals of different morphs, which usually coexist within populations at equal frequencies. Loss in the area and connectivity of suitable habitats may cause deviations from the expected equal morph frequencies. However, there is a need to evaluate the generality of this pattern at larger spatial extents and to identify possible underlying mechanisms. A citizen‐science approach was used to study morph frequencies in populations of the heterostylous grassland plant Primula veris across Estonia. We developed an online platform to facilitate an easy upload of the data. We examined the effect of the following variables in the surroundings of the study populations reflecting the landscape context on the deviation of morph ratios: (a) semi‐natural grasslands, (b) forests and shrubs, (c) human population density and (d) a proxy for plant population size. The citizen‐science approach provided unprecedented density of data from 1,700 localities. Nearly half of these observations, which were maintained for further analysis after data filtering, included over 62,000 short‐styled morphs and about 54,000 long‐styled morphs. Small populations were characterized by higher overall deviation of morph ratios from isoplethy (equal morph ratio). Deviation increased in semi‐natural grasslands located in regions with high human population density. The significant effect of human population density and plant population size on deviations of morph frequencies suggests the role of stochastic demographic effects of habitat fragmentation on morph ratios. Overall lower proportion of long‐styled morphs indicates that partial intra‐morph compatibility shown in long‐styled morphs may lead to higher inbreeding and related decline in fitness and abundance. Synthesis. Citizen‐science data about the morph type of Primula veris across Estonia obtained with the help of thousands of people demonstrates that in addition to plant population size, landscape context may affect plant reproductive traits, such as heterostyly. Larger population size of P. veris can help to buffer against random fluctuations in this trait. Increasing impact of human activities may have a negative impact on both small and large populations. The exact underlying mechanisms of the prevalence of one morph over the other, however, pose novel questions for further research.

Inter-morph differences are expressed also in other floral characteristics, such as the length of stigmatic papillae and the size of pollen grains. Charles Darwin in his landmark findings about the functional meaning of distyly in Primula veris was the first to suggest that such floral heteromorphism serves as an effective outcrossing mechanism (Darwin, 1862). Morphological dimorphism is accompanied by a genetically determined incompatibility system (Huu et al., 2016), which ensures crossing of plants with the ones characterized by the opposite style morph and prevents mating between individuals sharing the same floral morph characteristics. The heteromorphism is governed by a supergene; it comprises a series of genes closely spaced or linked on a chromosome. Short-styled individuals were thought to be heterozygous (Ss) and long-styled homozygous recessive (ss; Li, Webster, Furuya, & Gilmartin, 2007), but recent molecular genetic analyses in Primula suggest that the heterostyly-locus is a hemizygous region present in short-styled morphs, but absent from long-styled morphs (Huu et al., 2016;Nowak et al., 2015).
Rapid changes in land-use intensity, related negative effects on area and connectivity of natural and semi-natural habitats, and a loss of pollinator communities have brought the consequences of landscape change for heterostylous plants into the spotlight (Brys & Jacquemyn, 2015;Jacquemyn et al., 2012;Kéry, Matthies, & Schmid, 2003). The genetically determined incompatibility system in association with floral heteromorphism leads to symmetrical disassortative mating and generally equal morph frequencies at equilibrium conditions (Barrett & Shore, 2008). However, high spatio-temporal dynamics of natural and semi-natural habitats as a result of recent landuse change may shift such equilibria and related population viability in several ways. Severe habitat loss, being exemplified by the reduced area of European semi-natural grasslands, has substantially decreased the size of many plant populations (Kiviniemi, 2008). Resulting population bottlenecks and drift effects may lead to stochastic deviation of morphs from equal frequencies (Endels, Jacquemyn, Brys, & Hermy, 2002;Kéry et al., 2003). This, in turn, lowers the abundance of suitable partners for mating and thus reduces reproductive fitness Kéry et al., 2003). Nevertheless, the system of heterostyly is not completely leakage-proof (Brys & population size, landscape context may affect plant reproductive traits, such as heterostyly. Larger population size of P. veris can help to buffer against random fluctuations in this trait. Increasing impact of human activities may have a negative impact on both small and large populations. The exact underlying mechanisms of the prevalence of one morph over the other, however, pose novel questions for further research. Jacquemyn, 2015Jacquemyn, , 2020Wedderburn & Richards, 1990). According to Wedderburn and Richards (1990), crosses between L-morphs of P. veris may lead to seed set in up to 14.5% cases, while effective intramorph fertilization between S-morphs occurs very rarely despite the disproportionate inter-morph pollen flow in favour of short-styled morphs (Ornduff, 1980). Thus, L-morphs may have an advantage in conditions of the lack of compatible pollen (Van Rossum, De Sousa, & Triest, 2006) resulting from the scarcity of mates and/or an absence of pollinators, and should consequently be prevalent in fragmented habitats with smaller populations. With that, for example, high genetic diversity has been observed in populations dominated by L-morphs in distylous Pulmonaria officinalis, most probably owing to reproductive advantages of L-morphs in this species (Meeus, Honnay, Brys, & Jacquemyn, 2012). Yet, higher intra-morph fertilization between L-morphs can make this morph potentially more prone to biparental inbreeding (i.e. due to outcrossing between two L-morphs; Van Rossum & Triest, 2007), which may jeopardize the ability of populations to react to environmental changes.
Although previous evidence has demonstrated deviations of morph ratios from equal frequencies in small plant populations (e.g. Kéry et al., 2003), an understanding of the potential role of recent landscape dynamics and accompanying disruption of plant-pollinator interactions on morph ratios is limited (Barmentlo, Meirmans, Luijten, Triest, & Oostermeijer, 2018). While corresponding findings have attributed deviating morph ratios mainly to stochastic demographic processes, it is not clear whether landscape changes can induce directional deviations of morph frequencies as expected based on the differences in intra-morph fertilization patterns of S-and L-morphs.
In this study, we examine morph-ratio patterns in an extensive number of wild populations of the heterostylous grassland plant P. veris across Estonia. To obtain data from landscapes with different history of habitat fragmentation and to cover as wide a geographic scope as possible, we applied a citizen-science approach. In addition to enabling the collection of marked quantities of data, citizenscience projects are increasingly important in raising awareness of conservation issues (Ballard, Phillips, & Robinson, 2018;Turrini, Dörler, Richter, Heigl, & Bonn, 2018), an aspect crucial in an era of vast human-induced environmental change. Primula veris, whose abundance has shown a decrease due to recent land-use change (Brys & Jacquemyn, 2009), is thus a suitable species not only for gaining a better understanding of the phenomenon of heterostyly, but enables to increase knowledge of the importance of conserving semi-natural grassland habitats.
Over the last hundred years, the area of the primary habitat of P. veris (cowslip), i.e. semi-natural grasslands, has decreased more than 90% in Estonia with potentially detrimental effects on wild plant populations (Laasimer, 1965). We were interested in finding out whether such landscape changes affect the morph balance in the populations of P. veris. In particular, we hypothesized that the availability of semi-natural grasslands at the observed sites would yield more equal morph frequencies due to potentially larger population sizes and supporting pollinator communities. In contrast, declines in plant population size and higher proportion of unsuitable habitats replacing former grasslands, such as forests and regions with higher human population density, would lead to higher deviation of morph frequencies from equilibrium conditions. Finally, we hypothesized that in fragmented habitats, morph balance is skewed in favour of L-morphs due to partial intra-morph compatibility among L-morphs. F I G U R E 1 A screenshot view from the specifically designed web-platform (https://nurme nukk.ee/et) for an easy upload of the data within the frames of the citizen-science campaign 'Eesti otsib nurmenukke' carried out in Estonia in 2019. The view represents the explanation of the task of discriminating between short-styled S-morphs ('S-tüüp') and long-styled L-morphs ('L-tüüp') of Primula veris, and the tool to insert/remove data about the corresponding morph. Translation of Estonian text: 'Determine the morph of 100 plant individuals. In case there are fewer individuals, determine the morph of all individuals. In S-morph, several anthers are visible from above, while in L-morph only a single dot representing the stigma is visible' 2 | MATERIAL S AND ME THODS

| Study species
Primula veris L. (ssp. veris) is a perennial rosette hemicryptophyte typically occurring in calcareous semi-natural grasslands. Being a long-lived, insect-pollinated, habitat-specialist plant, P. veris falls into the category of species at a high risk due to factors of global change, such as habitat loss, related pollinator decline and climate change (Estrada et al., 2018). Previous studies have already demonstrated the negative impact of habitat fragmentation on the fitness and neutral genetic diversity of this species (Kéry, Matthies, & Spillmann, 2001;Van Rossum, Campos De Sousa, & Triest, 2004). However, owing to the relatively long life span of the species with an average length of 52.3 years (Ehrlén & Lehtilä, 2002), Lehtilä et al. (2016) observed significant time lags in the demographic patterns of P. veris in response to habitat deterioration. Similarly, P. veris showed considerable inertia to land-use change after grassland abandonment (Lindborg, Cousins, & Eriksson, 2005).
We considered P. veris to be a suitable species for the citizenscience campaign for several reasons. Primula veris is a relatively well-known species in Estonia, largely also because of a variety of uses of P. veris in folk medicine. Its peak of flowering arrives before most of the other plant species characteristic of semi-natural grasslands, which makes this species easily detectable in the landscape.
There are no other morphologically similar species flowering at the same time in Estonia and thus the risk of misidentification is low. The assignment of plants to L-and S-morphs is easy and can be reliably performed by anybody. However, in order to facilitate correct identification of the species as well as S-and L-morphs, the web-platform for uploading data was provided with detailed photos (Figure 1). Furthermore, participants were welcome to upload photos of their own observations, which enabled us to further verify the ability of observers to correctly fulfil the task. The species occurs all across Estonia, but is more widespread in northern and western parts of Estonia because of the higher proportion of suitable habitats (calcareous semi-natural grasslands) in this region.

| Data collection
The citizen-science campaign 'Eesti otsib nurmenukke' ('Estonia is looking for cowslips') aimed to engage participants from all over Estonia and a wide range of age groups. Because P. veris is flowering before summer holidays, we also invited schools and kindergartens to take part in the campaign. In collaboration with University of Tartu, Estonian Fund for Nature and design studio Fraktal, a userfriendly web-platform (www.nurme nukk.ee, now also available in English on platform www.cowsl ips.eu) was developed in Estonian and Russian to actively engage the Russian minority in Estonia (the native language of approximately 25%-30% of Estonian inhabitants is Russian). The platform enables easy uploading of all relevant data and GPS coordinates of the study locations with the help of smartphones. It was also possible to upload data via other devices with internet connection as well as to send data to the campaign organizers as a hard copy. The web-platform included detailed guidelines about carrying out the observations along with explanatory photos (Figure 1; Appendix S1). A more extended guideline in both Estonian and Russian was prepared for biology and kindergarten teachers to foster carrying out the task with a group of children (Appendices S2 and S3). In addition, the web-page included the fol- Each individual observation included the following baseline data: (1) Observation ID, (2) The date of observation, (3) Geographic location (the observer had a choice between automatic localization, writing the location text manually or assigning no location). In addition, the observers were asked the following questions: (4) Did you find cowslips? Option: Yes/No. (6) Observation of heterostyly. Participants were asked to report the morph type of a hundred individuals wherever it was possible, but were also welcome to submit data on a smaller/higher number of individuals depending on conditions (population size, possibility to spend time for a more extensive sampling effort).
Guidelines provided further details for carrying out the observations. In particular, participants were asked to determine the morph of individuals covering the whole area where cowslips were present (i.e. the whole population). After determining the morph of an individual, a participant had to go at least one step further (i.e. about 0.5-1 m) to avoid sampling the same individual. Web-assignment was provided with a clarifying photo and a short description of the particular morphological features and with '+' mark to report a corresponding morph ( Figure 1).
(7) Participants could upload up to five photos depicting the study site.
(8) People were welcome to add comments, e.g. a more detailed description of the study site, approximate population size and general feedback and suggestions.
(9) When interested in further feedback, participants could post personal details (the name of the observer, the name of the observing team, e-mail address for feedback correspondence).
(10) Final step concerned issues about privacy policy and the further use of personal data. It focused on the following key points: (a) the first citizen-science observation of the participant (agree/ do not agree), (b) additional data about nature can be sent to the before-mentioned contact information (agree/do not agree) and (c) the further use of privacy data (contact information, names) according to the conditions listed in the privacy policy and for receiving feedback about the outcome of the project (agree/do not agree). Uploading of heterostyly data did not depend on the agreement/disagreement of the participants with privacy policy.
However, we carefully considered these statements during further correspondence with participants and did not include any private data in analyses.
The web-platform and related guidelines were made public in the second half of April in 2019, i.e. the earliest possible date when wild populations of P. veris would start flowering in Estonia. Information about the citizen-science project together with guidelines for observations were advertised in a wide range of media platforms: national and local newspapers, Estonian Public Broadcasting TV news and other broadcasts, interviews in different radio programmes and social media (Facebook, Instagram) with live feed to the project web-page.
At the end of the campaign, the preliminary results of the project were introduced in outreach articles, in radio broadcasts and in feedback by e-mail to those participants who had agreed to this when filling protocol. In addition, to further acknowledge the contribution of participants, we organized a thanking event, where the most active participants and school teams were rewarded with various prizes.

| Data analyses
Plant population size has been shown to be a significant determinant of morph deviation from isoplethy (i.e. equal morph ratio; Kéry et al., 2003). Preliminary tests concerning the user-friendliness of the citizen-science platform suggested that the assessment of population size would involve too much uncertainty and potentially decrease the number of participants. Therefore, the task of evaluating population size was excluded from the protocol. However, in order to take into account the potential effect of population size, we used the number of observed individuals to create a categorical proxy variable with two levels indicating the abundance of plants at the observed location: (a) small populations (<100 individuals); (b) larger populations (≥100 observed individuals). The validity of this approach was further verified examining the comments, where many participants reporting the morph of fewer than required 100 individuals mentioned that there were <100 individuals at the observed site.
To examine the potential impact of habitat availability on morph ratios, we calculated the proportion of semi-natural grasslands in the landscapes surrounding those observations, which had been provided with GPS coordinates. The map layer of seminatural grasslands was obtained using the following sources: ( of semi-natural grasslands suitable for P. veris as well as the area of woody elements was assessed in circular buffers with radiuses of 100, 500, 1,000, 2,000 and 5,000 m to take into account both local and large-scale effects of landscape quality ( Figure 2).

F I G U R E 2
Location of the populations of Primula veris observed within the frames of the citizen-science project 'Eesti otsib nurmenukke' in Estonia, and an illustration for the buffer approach used for assessing the landscape characteristics (proportion of semi-natural grasslands and woody elements) in the surroundings of the study populations. The figure only includes those observations, which were provided with GPS coordinates  As we asked participants to consider carrying out the observa- were also interested in analysing the deviation from balance regardless of the directionality of the deviation (i.e. both towards S-or L-morphs), we also estimated the absolute value of SES S-morphs (SES ABS ) for each population.
In addition, to compare the main findings of this study to previous research, we assessed morph ratios using indices applied in previous studies of heterostyly in P. veris (i.e. Kéry et al., 2003). In particular, we assessed the proportion of S-morphs (DEV S-morphs ; note that Kéry et al. (2003) .
also estimated an averaged model based on Akaike weights, including the subset of models that were within 6 AICc points from the best-ranked model. To explore the effect of important predictors, we created partial regression plots for those predictors that had a significant (p < 0.05) effect on the averaged model. Finally, we estimated R 2 from this subset of models by using the Akaike-weighted averaged predicted values (Pistón et al., 2019). We checked spatial autocorrelation of the residuals of the averaged models by means of a Moran's I test. All statistical analyses were performed in r v3.6.0 (R Core Team, 2019).

| Absolute morph deviation
Human population density (cubic polynomial), the presence of semi-natural grasslands in the close surrounding of the populations (within 100 m) and their interaction were the most important determinants of the standardized absolute deviation from equal morph frequencies (SES ABS ) in populations of P. veris, and were included in all models with some support according to AICc (Tables 1a and 2).
These variables were the only ones with estimates that were statistically different from 0 in the averaged model. Population size (i.e. small and large) also had a relatively high Akaike weight, but its effect was not statistically significant. The rest of the predictors had lower values of importance with bedrock type and the proportion of grasslands within the radius of 5,000 m obtaining the lowest importance ( Table 2). The averaged (using Akaike weights) predictions of the selected models explained 17% of the variance in the dataset, and the residuals did not show significant spatial autocorrelation (Moran's I = −0.003, p = 0.819). SES ABS showed a marked nonlinear pattern with low values in the least densely populated areas and a noticeable increase for population densities higher than one person/ ha ( Figure 3).
The main determinants of DEV ABS were very similar to those of SES ABS with the main difference that population size played a more important role and was retained in all models selected according to AICc (Table 2). Smaller populations tended to deviate more strongly from even morph frequency than large ones ( Figure 4).

| Directional morph deviation
Population size, the presence of semi-natural grasslands at 100 m around the study populations of P. veris and human population density and its interactions with population size and presence of grasslands were the variables with the highest importance determining directional morph deviation SES S-morphs according to the model selection procedure (Table 2). None of the other predictors had substantial support with the presence of forests at 100 m around the populations being the next most important variable (importance = 0.52). Similar to the results of the models of absolute morph deviation SES ABS , bedrock type and the proportion of grasslands in the radius of 5,000 m gained the lowest importance in the models predicting SES S-morphs ( Table 2). The residuals of the averaged predictions did not show significant spatial autocorrelation (Moran's I = −0.007, p = 0.553). The proportion of S-morphs was higher in areas with high human population density, but only for populations occurring at semi-natural grasslands (Table 1b; Figure 5). Neither TA B L E 1 Averaged models (estimate, standard deviation and p-value) fitted for absolute (a; SES ABS and DEV ABS ) and directional (b; SES S-morphs and DEV S-morphs ) deviation of morph ratios in the study populations of Primula veris in Estonia. Variables: grassland (100 m)presence (1) or absence (0) of grassland areas in a 100 m buffer around the study population; proportion grassland (5,000 m)-proportional grassland area in a 5,000 m buffer (continuous value between 0 and 1); forest (100 m)-presence (1) or absence (0) of woody elements (forests and shrubs) in 100 m buffer; proportion forest (5,000 m)-proportional area occupied by woody elements in a 5,000 m buffer (continuous value between 0 and 1); bedrock (sandstone)-bedrock type at the observation location (takes value 1 for sandstone and 0 for limestone); small population: takes value 1 when the number of individuals sampled in the population is fewer than 100 and 0 otherwise; human density-human population density (continuous variable, number of persons per ha, log-transformed Directional morph deviation (DEV S-morphs ) and morph deviation accounting for the number of observed individuals (SES S-morphs ) were significantly correlated (r = 0.922, p < 0.001; Figure S2b).

| D ISCUSS I ON
In this study, we examined the effect of landscape context on the morph ratios of the distylous grassland plant P. veris in Estonia, where the area and connectivity of semi-natural grasslands has substantially decreased over the last century. A successful citizen-science campaign enabled us to obtain data from nearly 1,700 populations all across Estonia. S-morphs were more prevalent in the dataset compared to L-morphs. However, data from 740 observations, which met the criteria for further statistical analyses, revealed that landscape context had a significant effect on morph frequencies. In particular, morph TA B L E 2 Relative importance of the variables included in each model about the effect of human population density, population size (small/large) and landscape variables on the deviation of morph frequencies from isoplethy in populations of Primula veris in Estonia. For each explanatory variable we report the sum of Akaike weights of models including that variable, among the models within 6 AICc points of the best model, as well as the number of models that included it (in parenthesis)  (9) 0.09 (7) Bedrock:small population 0.05 (5) 0.04 (11) 0.14 (20) 0.17 (12) Forest (100 m

F I G U R E 3
The effect of human population density and its interaction with the presence of semi-natural grasslands at the study population or in the closest proximity (within the buffer with a radius of 100 m) on absolute morph frequencies in the populations of Primula veris in Estonia recorded in the frames of the citizen-science campaign. Summary of the fitted models is presented in Table 1a. R 2 represents the explained variance in the dataset according to the averaged (using Akaike weights) predictions of the selected models frequencies of P. veris were more likely to deviate in semi-natural grasslands located in regions with higher human population density ( Figure 3). Furthermore, the deviation in such grasslands was more often in favour of S-morphs ( Figure 5). We also found that smaller populations were more prone to overall morph deviations (Figure 4).

| Absolute morph deviation
We found that an increase in human population density in the surroundings of the observed populations caused a significant deviation of morphs from isoplethy ( Figure 3). In particular, plant populations in semi-natural grasslands located in regions with human population density higher than one person per hectare tended to have unbalanced morph frequencies. The average human population density for the whole of Estonia is relatively low with 0.31 persons/ha. Increasing human population density may indicate various processes in relation to anthropogenic pressure, such as urbanization and infrastructural development, but also higher rural human population in areas with high-intensity agriculture leading to the loss and fragmentation as well as decreased quality of habitats. Among other major shifts in landscape change in Estonia, such as the overgrowth of semi-natural grasslands with trees and bushes following the lack of management (Pärtel, Mändla, & Zobel, 1999), recent history has witnessed substantial urbanization and related urban sprawl in the surroundings of the larger settlements in Estonia. Suburban development has mainly occurred on former grasslands encompassing the suitable habitats of the study species P. veris (Roose, Kull, Gauk, & Tali, 2013).
Pressure from urban sprawl generally imposes a strong negative environmental impact, including the loss of natural and seminatural habitats, and lower quality of the remaining habitats (Hansen et al., 2005;Sushinsky, Rhodes, Possingham, Gill, & Fuller, 2013) leading to the loss of typical grassland species (Saar, Takkis, Pärtel, & Helm, 2012). Primula veris is a species characteristic to semi-natural grasslands, which in Estonia require continuous and regular moderate human disturbance through mowing and grazing for persistence. Indeed, moderate historic human population densities (most likely reflecting patterns of historic grassland management) was shown to have a positive effect on biodiversity (Helm et al., 2009;Pärtel, Helm, Reitalu, Liira, & Zobel, 2007).
However, in the same study system, contemporary human population densities had an adverse effect on the genetic diversity of a grassland-specialist plant (Helm et al., 2009). Similarly, our study indicates that contemporary human population density can skew patterns of morph deviations of P. veris in semi-natural grasslands.
Demographic instability as a result of abrupt plant population collapses (Kiviniemi, 2008) and related population bottlenecks bring along stochastic mortality of individuals, which can easily lead to deviated morph frequencies in populations having recently experienced loss in size. This has been confirmed by significant deviations of morph frequencies in small populations of P. veris (Brys, Jacquemyn, Endels, Hermy, & De Blust, 2003;Kéry et al., 2003), but also in other heterostylous species, such as P. vulgaris (Endels et al., 2002) and Pulmonaria officinalis (Brys, Jacquemyn, & Beeckman, 2008).

Models with different indices of deviation (SES ABS vs. DEV ABS )
led to relatively comparable results. However, overall deviation of morphs from balance (DEV ABS ; index sensu Kéry et al., 2003) depended also on population size: smaller populations (with fewer than 100 observed individuals) were characterized by higher deviation ( Figure 4). Interestingly, this finding was not confirmed for SES ABS , a measure standardized according to the number of observed individuals. Although frequencies of different morphs are much more prone to random fluctuation in smaller populations (Kéry et al., 2003), it is also probable that proportional deviation without standardization in relation to the number of observed individuals may not provide unbiased understanding of the mechanisms driving patterns of morph ratios. To take into account such biases, we mainly focused on standardized estimates of deviation (SES ABS and SES S-morphs ). R 2 = 0.14 Grassland 100 m absent Grassland 100 m present is why the overall dominance of S-morphs is thought-provoking.

| Directional morph deviation
Among previous research, only Lees (1971) recorded the dominance of S-morphs over L-morphs in P. veris and assigned this phenomenon to the potential advantage of heterozygous S-morphs in certain microclimates. Similarly, a sample of 522 individuals of P. veris collected by Charles Darwin exhibited an excess of 281 S-morphs over 241 L-morphs (Darwin, 1862). However, the number of sampled populations in those studies was too small to draw broader conclusions. In contrast, Kéry et al. (2003) who recorded the morph frequencies in 76 populations of P. veris in Switzerland, found no differences in the overall frequency of L-and S-morphs nor differences between the mean values per population.
One of the reasons for the overall prevalence of S-morphs could be the relatively higher fitness and survival of S-morphs owing to a stricter incompatibility mechanism among S-morphs and hence potentially higher genetic diversity. L-morphs, where a certain degree of within-morph compatibility has been observed (Wedderburn & Richards, 1990), may therefore exhibit higher (biparental) inbreeding, which can lead to lower survival of L-morphs. Indeed, Van Rossum and Triest (2007) reported higher kinship coefficients in L-morphs compared to S-morphs in fragmented populations of P.
veris lending support to our hypothesis about the relatively higher inbreeding in L-morphs. Preliminary genetic analysis of heterostyly data from fragmented semi-natural grasslands in western Estonia revealed moderate inbreeding in L-morphs, but surprisingly more so in populations dominated by S-morphs (Kaldra, 2020). It is, however, also possible that in such populations, high inbreeding in L-morphs has already led to the lower number of long-styled individuals compared to S-morphs. In the peripheral populations of P. merrilliana, Shao et al. (2019) observed intra-morph compatibility as well as self-compatibility not only in L-morphs, but also at similar frequencies in S-morphs. Furthermore, self-and intra-morph fertilization in S-morphs yielded a higher seed set than self-and intra-morph crosses in L-morphs of this species, which would help to buffer the increase of L-morphs as a result of intra-morph crosses (yielding to pins only). Thus, the variation of the patterns of intra-morph compatibility and self-compatibility in heterostylous plants may be higher than previously assumed and can also depend on the location of the population in relation to species' distributional range.
The drastic loss of semi-natural grassland habitats has disrupted vital plant-pollinator interactions necessary for successful reproduction (Valiente-Banuet et al., 2015). For outcrossing plant species depending on pollinator activity, the chronic lack of pollinators in combination with decreasing size of populations may favour phenotypes that are capable of partial self-compatibility or, in heterostylous species, partial intra-morph compatibility. Lower probability of legitimate pollen transfer can eventually lead to the loss of heterostyly (Barmentlo et al., 2018;Shao et al., 2019;Yuan et al., 2017). In our dataset, 53 observations (among those provided with GPS coordinates and with ≥10 individuals) included only one morph. These observations included fewer than 100 observed flowering individuals indicating small population size. Interestingly, 34 of those reported only L-morphs, while 19 observations comprised of only S-morphs, although the overall data included more S-morphs. Nearly all monomorphic observations originated from locations with no semi-natural grasslands nearby. Although monomorphism of P. veris in those locations needs to be further validated on site, these findings seem to suggest that the loss of suitable habitats may potentially indeed lead to drastic deviations of morphs with potential effects on fitness (Kéry et al., 2003) and evolutionary consequences, such as the loss of distyly. The dominance of long-styled monomorphism among our observations supports the hypothesis that partial intra-morph compatibility may to a certain extent be advantageous in alternative habitats of landscapes with no semi-natural grasslands, such as road verges, field boundaries and fallow land. The latter was also confirmed by the moderate prevalence of L-morphs in small plant populations in landscapes with higher human impact (Table 1b). However, what the longterm consequences of extreme deviations of morph frequencies are and whether some of these populations are indeed evolving towards the breakdown of heterostyly, needs further verification.

| Engaging people in citizen science
Citizen-science approaches for collecting environmental data have gained growing interest over the last decade (Silvertown, 2009;Turrini et al., 2018). Wider engagement of people in ecological research has multi-faceted benefits, from enhancing public knowledge in biodiversity and conservation issues to collecting scientific data at unprecedented spatial and temporal scales, which is often unattainable with traditional scientific methods (Ballard et al., 2018;Cooper, Dickinson, Phillips, & Bonney, 2007). Although we did not specifically assess the impact of our citizen-science approach on raising awareness of ecological knowledge and conservational issues, it was apparent from numerous comments uploaded by the participants that the impact was far beyond increasing knowledge just on heterostyly. In general, based on the comments, aspects covered by the citizen-science campaign apart from the main task of the campaign can be divided broadly into the following subtopics: landscape change, ecosystem services, botanical observations, interactions of observed P. veris individuals with other organism groups (herbivores and pollinators), habitats, where the observed populations were found (semi-natural grasslands, but quite often also in alternative habitats, such as road verges, quarries, former dump sites, airfield, graveyards etc) and environmental education. Thus, somewhat surprisingly, despite the rather specific task (i.e. determination of morph types in P. veris), people often approached the assignment more holistically. Comprehensive background information (e.g. explaining the wider implications of landscape change for biodiversity) provided in the web-platform as well as in media may have been among the potential reasons why people actually noticed more than asked to report in the data form.
Hence, such citizen-science projects may have wider implications than just enhanced knowledge on the particular subject.
Nearly 250 participating teams originated from schools or kindergartens. Positive feedback from teachers provided in the comments section of the protocol revealed that such participatory approaches offer very valuable educational experiences complementing the compulsory curriculum with both theoretical introduction into the subject as well as with related hands-on activities, which help to understand the topic. Furthermore, several teachers showed interest to participate in any follow-up projects on monitoring of heterostyly.
Understanding of the temporal dynamics of morph ratios is very scarce (Endels et al., 2002). Thus, both from educational angle and from the scientific point of view, continuation of the citizen-science project on heterostyly and expanding its geographic scope is surely worthy of further consideration.

| Implications for nature conservation
A number of P. veris populations observed in the current citizenscience campaign were located in regions exposed to relatively high human pressure, such as the county surrounding the largest town in Estonia, Tallinn. Urban sprawl and related human disturbances pose great challenges for maintaining biodiversity (Hansen et al., 2005). Morph deviation in the populations of P. veris appears to be to certain extent determined by stochastic demographic processes in fragmented habitats characterized by lower abundance of plant individuals and higher human population density. Furthermore, it is possible that the higher abundance of monomorphic populations with L-morphs reported in habitats with no semi-natural grasslands in the surroundings and the prevalence of L-morphs in small populations in landscapes with high human impact refers to the role of deterministic processes, such as the partial capability of L-morphs intra-morph fertilization. On the other hand, S-morphs were prevailing in most populations, particularly in semi-natural grasslands of landscapes experiencing higher human impact. Such patterns may affect the future viability of heterostylous species in various ways.
Previous research suggests that strong deviations in morph frequencies may further contribute to the loss of genetic variation (Meeus et al., 2012) and thus exacerbate the extinction vortex of heterostylous plants in fragmented habitats. Similarly, recent findings in western Estonia show that populations of P. veris with unbalanced morph frequencies exhibit lower genetic diversity (Kaldra, 2020).
Decreased genetic diversity in populations with deviated morph frequencies also make them more vulnerable to climate change and other environmental perturbations (Leigh, Hendry, Vázquez-Domínguez, & Friesen, 2019).
Novel selection pressures, such as the scarcity of suitable mates and lack of pollinating insects, may eventually lead to both stochastic and directional deviations of morphs and, ultimately, to evolutionary changes in reproductive traits . In addition to the revealed effects of landscape context on morph deviations in P. veris, this study yielded a number of observations with only a single morph present, particularly in small populations.
These findings need to be examined further in order to clarify what the mechanisms for the unbalanced morph frequencies are and whether the development of homostyly as a result of recent landscape change is more common than observed until now.

ACK N OWLED G EM ENTS
We are very grateful to all people who provided data on heterostyly within the frames of the citizen-science campaign 'Eesti otsib nurmenukke' in 2019. We thank people in Estonian Fund for Nature, and also others who contributed to the organization of the campaign, website designers from design studio Fraktal for creating an eyecatching and easy-to-use web-platform and science communicators for delivering information about the campaign. We are particularly grateful to teachers who engaged numerous pupils in observations of heterostyly and provided valuable feedback. We thank two anonymous reviewers and the Associate Editor for very useful and inspiring comments about the first version of the manuscript. We

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/1365-2745.13488.