Habitat segregation of plate phenotypes in a rapidly expanding population of three-spined stickleback

. Declines of large predatory ﬁ sh due to overexploitation are restructuring food webs across the globe. It is now becoming evident that restoring these altered food webs requires addressing not only ecological processes, but evolutionary ones as well, because human-induced rapid evolution may in turn affect ecological dynamics. We studied the potential for niche differentiation between different plate armor phenotypes in a rapidly expanding population of a small prey ﬁ sh, the three-spined stickleback ( Gasterosteus aculeatus ). In the central Baltic Sea, three-spined stickleback abundance has increased dramatically during the past decades. The increase in this typical mesopredator has restructured near-shore food webs, increased ﬁ lamentous algal blooms, and threatens coastal biodiversity. Time-series data covering 22 years show that the increase coincides with a decline in the number of juvenile perch ( Perca ﬂ uviatilis ), the most abundant predator of stickleback along the coast. We investigated the distribution of different stickleback plate armor phenotypes depending on latitude, environmental conditions, predator and prey abundances, nutrients, and benthic production; and described the stomach content of the stickleback phenotypes using metabarcoding. We found two distinct lateral armor plate phenotypes of stickleback, incompletely and completely plated. The proportion of incompletely plated individuals increased with increasing benthic production and decreasing abundances of adult perch. Metabarcoding showed that the stomach content of the completely plated individuals more often contained invertebrate herbivores (amphipods) than the incompletely plated ones. Since armor plates are defense structures favored by natural selection in the presence of ﬁ sh predators, the phenotype distribution suggests that a novel low-predation regime favors stickleback with less armor. Our results suggest that morphological differentiation of the three-spined stickleback has the potential to affect food web dynamics and in ﬂ uence the persistence and resilience of the stickleback take-over in the Baltic Sea. and topographic openness, in bays dominated by stickleback (a) or perch (b). The solid trend line shows a signi ﬁ cant relationships between the number of stickleback and topographic openness in the bays dominated by stickleback, as determined by linear regression (df = 14, r 2 = 0.37, t = 2.9, P = 0.013).


INTRODUCTION
Human-induced biodiversity loss strikes particularly hard on higher trophic levels, and many populations of top predators have declined across the globe (Ripple et al. 2014). In marine systems, loss of larger predatory fish has induced changes in predator-prey relations that propagate through food webs, affecting biodiversity, and change ecosystem structure (Jackson et al. 2001, Daskalov 2002, Worm et al. 2006. In response to the severe economic and cultural losses resulting from these changes, various measures have been initiated to restore biodiversity at higher trophic levels (e.g., marine protected areas; Gell andRoberts 2003, Lester et al. 2009). However, these measures may be ineffective because of pervasive but poorly understood feedback mechanisms within the ecosystem (Nystr€ om et al. 2012), as well as interactions between multiple human impacts that delay or prevent recovery (Frank et al. 2005, 2011.
Currently, we have an inadequate understanding of the drivers that structure food webs that are, themselves, changing. One poorly understood factor is rapid evolution. Evolution was traditionally assumed to be slow and unimportant on ecological timescales and has therefore rarely been included in ecological studies or ecosystem management (Fussman et al. 2007, De Meester et al. 2019, Hendry 2019. Today, however, we know that selection can change population dynamics on short timescales and that overfishing in particular has the potential to cause exceptionally rapid evolution of target species (Conover and Munch 2002, Kuparinen and Meril€ a 2007, Hutchings and Fraser 2008, Audzijonyte et al. 2013. Over the past decade, we have also become increasingly aware that rapid evolution can feed back on to ecosystem properties, so-called eco-evolutionary feedbacks (Post and Palkovacs 2009, Schoener 2011, Turcotte et al. 2011, Hendry 2017, Govaert et al. 2019. Still, evolutionary consequences of large-scale predator loss for lower trophic levels are not widely understood.
Declines of large predatory fish are typically followed by dramatic increases in the abundance of their main prey, often smaller fish species; a process known as mesopredator release (Prugh et al. 2009). For the prey fish, the strong increase in density intensifies intraspecific competition and may shift the selection regime from predation defense toward increased individual-level resource specialization and population-level resource use diversity (Svanb€ ack and Persson 2004, Svanb€ ack and Bolnick 2005. Accordingly, there are clear examples of how populations have increased their competitive ability and diversified after escaping control from their native predators (Palkovacs et al. 2011, Mlynarek et al. 2017. Thus, mesopredator release caused by the loss of large predators may lead to a strong selective pressure for different traits among prey fish. This may, in turn, impact other trophic levels by introducing new predator-prey interactions, ultimately impacting the nature and strength of trophic cascades. Here, we explore the potential for niche differentiation in a population of prey fish in the Baltic Sea, which have increased dramatically the past decades at the same times as a number of predator populations have decreased strongly.
In the western parts of the central Baltic Sea (NW Europe), densities of a small mesopredatory fish-the three-spined stickleback (Gasterosteus aculeatus L., hereafter stickleback)-have increased up to 45 times during the last decades, following regional and/or local declines of their natural predators; the large-bodied commercially and recreationally harvested fish cod (Gadus morhua), pike (Esox lucius), and perch (Perca fluviatilis; Ljunggren et al. 2010, Eriksson et al. 2011, Bergstr€ om et al. 2015, Olsson 2019. Suggested mechanisms behind the stickleback population explosion include a combination of human impacts that link drivers of change across large spatial scales. Juvenile stickleback migrate offshore during autumn and, when mature, return to the coast for spawning in spring (Bergstr€ om et al. 2015). Offshore, overfishing of cod, increased water temperatures, and increased food availability due to eutrophication may have contributed to their increased survival (Eriksson et al. 2011, Lefebure et al. 2014. Inshore, local decreases in perch and pike populations may have increased recruitment and juvenile survival of stickleback in some areas ( € Ostman et al. 2014, Bergstr€ om et al. 2015, Hansson et al. 2018. The local decreases in pike and perch are caused by loss or degradation of reproduction habitats due to wetland drainage, coastal constructions, and boating (Nilsson et al. 2014, Sundblad and Bergstrom 2014, Hansen et al. 2019) and, potentially, increased predation from seals and cormorants ( € Ostman et al. 2014, Hansson et al. 2018, Veneranta et al. 2020. Today, there is a strong inverse relationship between local abundances of stickleback on the one hand and perch and pike on the other (Bergstr€ om et al. 2015), and we have indications that in large abundances the stickleback themselves limit the perch population by consuming their eggs and larvae (Bystr€ om et al. 2015, Ekl€ of et al. 2020. Spatially explicit analyses demonstrate that a shift in dominance from perch and pike to stickleback-dominated bays has propagated from the outer archipelago toward the coast as a spatially propagating wave through time (Ekl€ of et al. 2020). The dramatic increase in stickleback abundance, in combination with elevated resource availability from large-scale eutrophication, has restructured near-shore food webs and increased filamentous algal blooms; thereby, threatening coastal biodiversity (Eriksson et al. 2009, Sieben et al. 2011a, Bystr€ om et al. 2015, Donadi et al. 2017. Globally, the marine stickleback is known to have rapidly and repeatedly adapted to freshwater habitats with low-predation pressure by evolving reduced armor phenotypes (i.e., lower number of bony lateral plates and loss/reduction of pelvic girdle), apparently in response to relaxed selection on these antipredator traits (Bell 2001, Huntingford andCoyle 2006). A large body of research supports the hypothesis that fully plated marine individuals are better protected against predation by fish than individuals with reduced plating (e.g., Reimchen 1992Reimchen , 2000. On the other hand, reduced lateral plate numbers can increase fast-start performance and maneuverability and thereby provide means to avoid predation in habitats where refuge from predators (e.g., vegetation) is present (Bergstr€ om 2002, Leinonen et al. 2011b). In addition, adaptation to predation can strongly influence the feeding ecology of stickleback (reviewed in Bell andFoster 1994, Harmon et al. 2009).
Three lateral plate phenotypes are described for the three-spined stickleback in the Baltic Sea: fully plated individuals, which have plates along the sides of the body that run uninterrupted from the head to the tail fin (caudal end); partially plated individuals, which lack plates on the midsection of the body; and low plated individuals, which only have plates on the head (pectoral) region of the body (M€ unzing 1963(M€ unzing , Reimchen 2000. The plate phenotype is highly heritable and largely controlled by one major gene (Ectodysplasin; Colosimo et al. 2004), and hence, indicative of the underlying genotype (Loehr et al. 2012). In the Baltic Sea, all three plate phenotypes of the three-spined sticklebacks have co-occurred for at least a century (Aneer 1974, and there is evidence for adaptive phenotypic differentiation in number of lateral plates across the Baltic Sea . The aim of this study was to evaluate the potential for niche differentiation between different plate morphs as a response to changing predation regimes, in the dramatically expanding population of three-spined stickleback in the Baltic Sea. We first assessed temporal trends in the abundance of the three-spined stickleback and juvenile perch, the most common predator, using time-series data covering 22 years . Then, we related differences in the spatial distribution of plate phenotypes to the abundance of predatory fish and production of the community. Finally, we documented differences in the diets of the different plate phenotypes, to evaluate the potential for feedbacks on food web structure caused by changes in plate phenotype frequencies. We hypothesized (1) that in areas with relatively low abundances of predatory fish, reduced selection for antipredator traits results in an increased proportion of incompletely (low or partially) plated stickleback, compared with completely plated individuals; and (2) that the v www.esajournals.org different stickleback plate phenotypes have different feeding ecologies, as reflected by different prey compositions in their stomach contents.

MATERIAL AND METHODS
Trends in three-spined stickleback and juvenile perch abundances over time We first analyzed temporal trends in the abundance of juvenile three-spined stickleback and perch (the dominant coastal predator on stickleback). Data on temporal trends of juvenile threespined stickleback and perch abundances from 1996 to 2017 were extracted from the Swedish national database for coastal fish (http://www. slu.se/kul). The data combined samplings of juvenile fish from various monitoring programs and research projects that quantified fish recruitment along the Swedish Baltic Sea coast. The available data on juvenile fish are spatially far more comprehensive than those for adults. Therefore, use of juvenile data enabled us to extract a time series that covered the full outer to inner archipelago gradient, which better represents the full range of habitats occupied by different species of fish than the regular coastal fish monitoring programs (see, e.g., HELCOM 2018). The database on juvenile fish also includes small species, such as stickleback, which are poorly represented in the regular monitoring programs (as the multimesh gillnets used for monitoring lack mesh sizes small enough to capture stickleback in a representative way; HELCOM 2018).
We extracted abundances of juvenile perch and three-spined stickleback from an area covering a 360 km stretch of the central Baltic Sea coastline (counties of Uppsala, Stockholm, S€ odermanland and € Osterg€ otaland; from Skutsk€ ar to V€ astervik). This area overlapped with a separate sampling of stickleback phenotypes in 2014, presented below. The data on juvenile fish collected in the Swedish national database for coastal fish, were sampled in shallow bays using small underwater detonations; a standard method that catches fish with a swim bladder of lengths up to 20 cm within a blast radius of 2.5-5.0 m depending on detonation strength (1-10 g; size range of fish 2-10 cm; Snickars et al. 2007). Since fish are heterogeneously distributed at small spatial scales and the detonation method takes a snapshot of a defined area (~20-80 m 2 depending on detonation strength), the abundance data is highly variable. In addition, bays have been sampled at various intensities (from 5 to 90 detonations per bay, mean AE SD = 16 AE 11). We therefore square-root transformed abundance data from each detonation and then used averages for each bay for the statistical analyses. Moreover, we only included the years 1996-2017 in the analysis because before 1996 most years were only represented by one single bay; leaving 7553 detonations from 480 bayyear combinations in the dataset. The mean depth of all sampling was 1.7 AE 0.59 m (AESD) and nearly all sampling (94%) was conducted during late July-early September. All sampling was performed by certified personnel and with required ethical permits and exemptions from national fishing regulations.
We analyzed trends in stickleback and perch abundances by fitting linear models with year as an explanatory variable. We did not perform a formal time-series analysis accounting for autocorrelation of sampling stations over time, because most of the bays were sampled only once (184 out of 256 unique bays) and less than 10% were sampled more than four times. In addition, the replication of bays for each year was highly skewed: The first five years had an average replication of~10 bays per year, compared with~30 for the following decades. To balance the influence of the different time-periods on the model fit, we therefore used the averages for each year as response variable and treated each year average as a true replicate. However, this reduces the variability of the data, meaning that the R 2 of the model fit only is related to the variation in yearly means; not to the variation between bays within or across years. Since the sub-set of bays were different each year, we included a spatial structure in the model to make sure trends were not caused by differences in sampled regions or characteristics of bays. The spatial structure consisted of adding average sampling depth, latitude, and longitude as explanatory factors to the model. We selected the best model for perch and stickleback abundances based on AIC c (AIC for small samples); using the MuMIn package in R 3.5.2 (version 1.43.15; R Core Team 2018). Two years (1996 and 2000) were excluded from the model selection because they caused strong collinearity between latitude and longitude; sampling predominantly occurred in the southern and innermost part of archipelago.
v www.esajournals.org After the selection procedure, the omitted years were reintroduced into the final model. In the reduced dataset there was no correlation between the explanatory variables and we did not detect any problems with collinearity in the candidate models (the variance inflation factor [VIF]: year = 1.10, latitude = 1.19, longitude = 1.22, depth = 1.06), or the final models (VIF: year = 1.00, depth = 1.00). If necessary, variables were normalized using square-root transformations.

Characterizing the coastal food web and stickleback phenotype composition
Since the time-series data described above mainly contain juvenile fish abundances, we tested which factors that best described stickleback abundances and plate phenotype distribution by characterizing the benthic food web in 32 shallow coastal bays in spring 2014 (max sampling depth between 3 and 6 m). The area covered was the same as where the time-series data was extracted from (Appendix S1: Table S1, Fig. S1). The bays represented a gradient from inner to outer archipelago, characterized by gradients in wave exposure and bay morphology; factors that regulate the composition of the biological communities in the Baltic Sea archipelago (see Appendix S1 for detailed explanation of selection criteria and calculations). Some of the field data were included in a structural equation model published in Donadi et al. (2017), which specifically analyzed top-down effects of sticklebacks on lower trophic levels. Here we focus on explaining the distribution of the different plate phenotypes of stickleback-which is previously un-published data.
The fish community was sampled in spring (from May to early June) using three to six Nordic survey gillnets set out overnight in each bay (mesh size 5-55 mm, details in supplementary information in Appendix S1). All fish caught were identified to species, measured for total length, and counted. Bleak (Alburnus alburnus), ruffe (Gymnocephalus cernuus), three-spined stickleback, perch, and roach (Rutilus rutilus) dominated the fish fauna (Appendix S1: Table S2). Bleak, perch, ruffe, roach, and three-spined stickleback all eat a wide variety of food items, including plant material, invertebrates, and zooplankton; but larger perch also more commonly include medium-sized fish in its diet (Appendix S1: Table S2; Jacobson et al. 2019). The number of individuals of each fish species caught was pooled for each bay and catch per unit effort (CPUE) was expressed as number of individuals per net and night of each species. In each bay, we estimated vegetation and filamentous algae cover, collected zooplankton, and measured depth, fluorescence, temperature, and turbidity, at three to six randomly distributed stations (Appendix S1: Fig. S2a). Mesograzer (benthic macro-herbivores) biomass was quantified at each station by sampling epifauna in a randomly placed 0.2 9 0.2 m frame connected to a 1 mm-mesh bag (Appendix S1: Fig. S2b,c). After identification of mesograzers in the laboratory (i.e., shredders and gatherers, whose diet is dominated by macrophytes: see Donadi et al. 2017), the biomass of mesograzers was estimated as gram ash-free dry-mass (AFDM) using taxonspecific correlations between length and weight (Ekl€ of et al. 2017). The number of both stations and nets depended on the size of the bay. They were randomized within depth strata from a topological map to be representative of the depth profile in each bay, but restricted to a depth between 0.5 and 3 m (Appendix S1; also see Donadi et al. 2017 for a detailed description of sampling procedures).
In each bay, we sub-sampled 30 individuals of the sticklebacks caught in the nets for morphological and diet analyses. However, in 16 of the bays we caught less than 30 individuals (mean = 5.5, SD = 5.3, min = 0, max = 22; CPUE). In the other 16 bays we caught high abundances of stickleback and the fish community was often completely dominated by three-spine stickleback individuals (stickleback bays: mean = 860.3, SD = 662.4, min = 32, max = 1799; CPUE). Earlier studies show an inverse relationship between stickleback and perch abundances (Bergstr€ om et al. 2015, Ekl€ of et al. 2020. Accordingly, the average abundance of perch was three times higher in the 16 bays where we found less than 30 stickleback (mean = 87.1, SD = 65.0; CPUE), in comparison to the 16 bays where we caught more than 30 stickleback (mean = 28.3, SD = 21.6; CPUE). To classify the bays, we calculated a mesopredator index as the total number of sticklebacks in each bay divided by the total number of perch and stickleback together. Thus, if the bay had stickleback present but no perch, v www.esajournals.org the index equaled 1; if the bay had perch present but no stickleback, the index equaled 0. The index confirmed the negative relationship between perch and stickleback (Appendix S1: Fig. S3): The 16 bays where we found less than 30 stickleback had an index of less than 0.4 (perch-dominated bays) and the 16 bays where we caught more than 30 stickleback had an index of more than 0.6 (stickleback-dominated). To get a reliable and comparable estimate of the fraction of plate phenotypes, we only included the 16 sticklebackdominated bays where we could collect ≥30 individuals in the analyses of plate phenotype distributions (Appendix S1: Table S1). To determine body shape variation, eleven length parameters were documented from all individuals in the stickleback-dominated bays (following Jones et al. 2012, Appendix S2). The number of lateral plates was counted on both the left and right sides of the body, and the average number of lateral plates for each individual was calculated. Plates were counted starting immediately after the operculum to the end of the caudal peduncle, following Aneer (1974). Plate phenotypes were identified by looking at gaps in the plate armor: Low plated fish lacked caudal plates; partially plated fish had a clear gap between the main body plates and the caudal plates (the width of at least two body plates); fully plated fish had no clear gap in the plating. However, there were a number of phenotypes in the transition phase between low and partially plated individuals; for example, individuals with only a few clear body plates followed by a clear gap and then a keeled tail consisting of reduced plates. For the analyses, we therefore divided the phenotypes into completely (fully) or incompletely (low and partially) plated sticklebacks. In total, we morphotyped 492 individuals from the 16 stickleback-dominated bays. Since we only defined the stickleback as completely or incompletely plated, their relative proportions of the phenotypes are the inverse of each other. In the analyses, we therefore only modeled the proportion of incompletely plated three-spined stickleback individuals.
Algal production was estimated for each bay as percent cover of filamentous algae on settling plates following three months in the water (May-August). Each plate consisted of four 5 9 5 cm unglazed ceramic tiles glued in pairs to two bricks, which were placed flat on the bottom at~1.5 m depth (Appendix S1: Fig. S2d). Algal cover on settling plates is a good measure of relative net primary production in the Baltic Sea (e.g., Eriksson et al. 2006Eriksson et al. , 2009Eriksson et al. , 2012, as filamentous algae respond to nutrient enrichment by fast accumulation of biomass, but are simultaneously grazed by invertebrate mesograzers (Worm et al. 2002, R aberg andKautsky 2007). Mesograzers are in turn eaten by fish, and thereby play a key role in trophic cascades (Eriksson et al. 2009, Donadi et al. 2017. In eight of the 32 bays, individual tiles were either lost or completely/partially covered with sediment. These bays were therefore not included in the analysis of algae cover development. We characterized the physical environment of each bay by estimating depth, wave exposure, and topographic openness (Ea). Depth was estimated as the average depth of all sampling points in each bay. Wave exposure was estimated for each bay using a simplified wave model based on fetch, topography, and wind conditions using Geographical Information System (GIS) methods ). Ea is defined as 100 9 At/a: where At is the crosssectional area of the smallest section of the opening of the bay towards the sea and a is the bay surface area. Topographic openness is one of the most important factors structuring the coastal communities of Baltic Sea bays (Hansen et al. 2012). We used topographic openness as a proxy for production in the bays, because we wanted a factor in the model that is independent of stickleback dynamics. The more direct measures, the development of algae and vegetation, are influenced by indirect top-down and direct fertilization effects of the stickleback themselves (Sieben et al. 2011a,b, Donadi et al. 2017. With increasing openness the influence of seawater and the influx of nutrients from surrounding water bodies increases. This promotes the benthic community (Donadi et al. 2017), which is reflected in higher biomass of marine and brackish water algae and crustacean macroherbivores (mainly benthic amphipods) in more open bays (Hansen et al. 2012). We tested our assumption that topographic openness determine production by modeling the relation with vegetation cover and algae development on the settlement tiles (linear models; lm function in R, version 3.6.1).
We tested which factors that best described stickleback abundances (# CPUE) and plate phenotype distribution (proportion of partially plated) by comparing General Linear Models including explanatory factors related to spatial structure (latitude, longitude), physical characteristics (wave exposure; depth [average depth of sampling points in bay]; and topographic openness [Ea]), prey availability (mesograzers and zooplankton), water temperature (at sampling), nutrient availability (the limiting nutrient phosphorus), and predation pressure (cpue perch >25 cm). We used perch abundance as an indicator of predation pressure because perch is by far the most common stickleback predator in the coastal fish community. We only included perch >25 cm, as stickleback then constitute their most common prey item (>40% of the diet in spring; Jacobson et al. 2019). The other common coastal fish predator, pike, is poorly represented by the gill-net catches because of their sedentary lifestyle, and was therefore not included in the analyses. We selected the best linear model (lm function) for stickleback abundances and plate distribution based on AIC c (AIC for small samples); using the MuMIn package in R 3.5.2 (version 1.43.15; R Core Team 2018). Data were square root ( stickleback; Ea; # perch >25 cm), log10 (# zooplankton; # mesograzers; temperature) or logit (proportion partially plated stickleback) transformed, to fit assumptions of normality. The logit transformation of plate phenotype was calculated as the natural logarithm of the proportion of partially plated individuals divided by the proportion of fully plated individuals (LN[proportion partially plated/proportion fully plated]).

Stickleback diet analyses
We used a metabarcoding approach to analyze diet differentiation between sympatric stickleback plate phenotypes. In many of the bays with few sticklebacks, we only found one of the two plate phenotypes. To avoid bias in stomach content analyses caused by the food availability in those bays with only one phenotype present, we randomly selected 177 individuals from 10 bays where both plate phenotypes were present (Appendix S1: Table S1). In the laboratory, the stomachs were dissected and flushed with 80% EtOH to remove all stomach contents and stored at À20°C in 80% EtOH. DNA was then extracted from the stomach contents using the UltraClean Tissue and Cells DNA Isolation Kit (MO BIO Laboratories). The dual PCR amplification method was then used for Illumina MiSeq library preparation (Bourlat et al. 2016). The amplicon primers were based on Leray et al. (2013) yielding a 313 bp fragment targeting the cytochrome c oxidase subunit 1 mitochondrial gene (CO1), and a blocking primer for G. aculeatus was also used to prevent amplification of predator DNA (all primer sequences can be found in Jakubavi ci ut _ e et al. (2017). Illumina MiSeq produced 30,103,790 paired-end reads of 300 bp in length. The processing steps were performed using Qiime (Quantitative Insights into Microbial Ecology) version 1.9.1 (Caporaso et al. 2010) and custom python scripts. Paired-end joining was done using the Qiime fastq-join tool. Primer sequences were removed using a custom python script, and remaining chimeric reads were excluded using UCHIME (Edgar et al. 2011). The Bayesian clustering algorithm CROP was used to cluster the sequences into operational taxonomic units (OTUs; Hao et al. 2011). Taxonomic assignment was performed with the Uclust software implemented in Qiime (Edgar 2010), using a 97% similarity limit when comparing the CO1 sequences with our own reference database of Chironomidae, Nemertea, Xenacoelomorpha, and Oligochaeta, combined with barcodes of Swedish Echinodermata, Mollusca, Cnidaria and Arthropoda from the Swedish Barcode of Life database (SweBol). Detailed methods for DNA extraction, amplicon library preparation, and bioinformatic analyses of the current samples can be found in Jakubavi ci ut _ e et al. (2017). The raw sequence data are available as fastq files from the NCBI sequence read archive under accession number SRP101702, BioProject number PRJNA378633.
For analysis, the OTU tables were converted to a presence/absence matrix. The prey species were divided into broad taxonomic groups relevant for the ecology of the system: crustacean mesograzers (amphipods, isopods, and mysids), mollusks (gastropods and bivalves), benthic worms (polychaetes and annelids), insects (chironomids and others), and zooplankton/meiofauna (cladocerans, copepods, and ostracods). For the complete list of prey taxa found in this v www.esajournals.org study, see Table S1 in Jakubavi ci ut _ e et al. (2017). A table of presence or absence of prey groups in the stomach of each fish was compiled. Differences in the occurrence of prey in the stomachs of completely compared to incompletely plated phenotypes, was analyzed using the multivariate statistics package Vegan (version 2.5-6) in R. First, we performed an Analysis of Similarities (ANOSIM) to test if there was a difference in the composition of prey in the stomachs of completely and incompletely plated sticklebacks (model: species composition of prey~plate phenotype). Then we plotted prey groups and the centroids of plate phenotypes using the information from an unconstrained NMDS ordination, to identify which prey groups contributed most to the differences between plate phenotypes. For each of the prey groups that was identified by the visual inspection (amphipods, mysids, bivalves, and polychaetes; see Results) we analyzed the difference in stomach content between completely and incompletely plated fish, using a generalized mixed-effects model using the lme4 package in R (version 3.5.2; Bates et al. 2015). The model was based on the binomial error distribution with a logit link function (logistic regression), and included plate phenotype as a categorical explanatory variable, and bay as random factor.

Changes in predator-prey relations along the coast
Juvenile three-spined stickleback abundances increased strongly along the western coast of the central Baltic Sea between 1996 and 2017 (Fig. 1). The best fitting model only included year as an explanatory factor (exponential fit: F 1,20 = 42.2, R 2 adj = 0.66, P < 0.001). In contrast, juveniles of the dominating coastal predator perch decreased linearly between 1996 and 2017 (linear fit 1996-2017: t = À2.3, P = 0.034), and sampling depth was included as a very important factor in the best fitting model (abundance increased with sampling depth: t = 4.1, P < 0.001; full model: F 2,17 = 12.2, R 2 adj = 0.54, P < 0.001; Fig. 1). Thus, the relative abundances of juveniles of top-and mesopredators have changed dramatically over the past decades, likely reflecting a changed trophic structure along the coast.

Stickleback abundances and the distribution of plate phenotypes along the coast
The food web field survey indicated the importance of production and predation in regulating the three-spined stickleback population. Stickleback abundance was best described by a model including topographic openness and the number of large predatory perch (model recommended by AIC c : R 2 adj = 0.43, F 1,29 = 12.7, P < 0.001; Appendix S3: Table S1), where stickleback abundance increased with increasing topographic openness (Ea: df = 29, t = 4.2, P < 0.001) and decreased with increasing abundance of large predatory fish (number of perch >25 cm: df = 29, t = À1.9, P = 0.065). There was no correlation between the explanatory variables topographic openness (Ea) and CPUE of large perch (Pearson's product moment correlation: r = À0.19, t = À1.03, df = 30, P = 0.309), and we did not detect any problems with collinearity in the model (VIF = 1.04). However, the individual effect of perch was only a statistical trend and the fit of the residual after correcting for topographic openness was poor (residual stickleback [sqrt CPUE~Ea] depending on # perch >25 cm: R 2 adj = 0.08). Considering the bimodal distribution of bays, where half where dominated by perch and the other half by stickleback, we further explored the data by allowing for an interaction between topographic openness and bay type (perch or stickleback dominated) instead of perch numbers. This model improved the fit considerably and showed that stickleback increased strongly with topographic openness in stickleback-dominated bays, but not in perchdominated bays (full model: R 2 adj = 0.77, F 3,28 = 12.7, P < 0.001; significant interaction effect: df = 28, t = 2.3, P = 0.028; Fig. 2). Topographic openness significantly explained estimates of net production in the bays; both general vegetation cover (R 2 adj = 0.22, F 1,30 = 9.6, P = 0.004, Fig. 3a) and algal cover development on the settlement tiles (R 2 adj = 0.25, F 1,23 = 9.1, P = 0.006; Fig. 3b) increased with openness.
We observed two distinct stickleback phenotypes corresponding to completely (fully) and incompletely (partially/low) plated phenotypes and very few intermediates (Fig. 4). The plate phenotype distribution was best described by a model including topographic openness and the number of large predatory perch (number of perch >25 cm; model recommended by AIC c : R adj = 0.39, F 2,13 = 5.9, P = 0.015; Fig. 5; Appendix S3: Table S2), where the proportion of incompletely plated individuals increased with topographic openness (t = 2.7, n = 16, P = 0.019) and decreased with the abundance of large perch (t = À2.5, n = 16, P = 0.026). This indicates a potential difference in selective regime on plate phenotypes depending on bay characteristics and predation. The incompletely plated phenotype was five times more common than the completely plated phenotype, and their relative abundance increased from 75% in the most enclosed bays to 90% in most open bays (Fig. 5a). At the same time, bays where perch was more abundant had a lower proportion of incompletely plated stickleback than expected from the topographic openness of the bay (Fig. 5b).

DISCUSSION
Our results suggest that decreasing predator populations influence the selection regime of three-spined stickleback across the western Baltic Sea. Coastal time-series data show that a dramatic increase in stickleback abundance coincided with a decline in the abundance of juveniles of its main predator, perch, between 1995 and 2015. Suggested causes for the decline of perch include loss of recruitment habitats (Sundblad and Bergstrom 2014), but recent studies indicate that sticklebacks can also control the recruitment of perch by feeding on its juvenile stages (Bystr€ om et  increased with topographic openness and decreased with the abundance of large predatory perch. The relation between the distribution of plate armor phenotypes and predatory perch indicates that declining predator populations have the potential to enable a shift in selection pressures for the increasing stickleback population: from high predation pressure to high intraspecific competition. Moreover, the two plate armor phenotypes had small but significant differences in their diets. This indicates that niche differentiation between the plate phenotypes may even alter the effects of the stickleback population on food web dynamics, but this is a hypothesis that needs to be further explored. The frequency of incompletely plated stickleback phenotypes increased with increasing topographic openness and decreasing abundances of large predators. The relation with topographic openness may depend on the general increase in production with openness in this system (Hansen et al. 2012, Donadi et al. 2017; meaning increasing abundances of invertebrate and plankton prey reducing intraspecific competition, but also increasing vegetation cover (see Fig. 4a). Vegetation cover may favor fewer armor plates since unarmored sticklebacks may move through vegetation cover more effectively (Bergstr€ om 2002, Leinonen et al. 2011b. The relation between plate phenotype and predation is interesting because the three-spined stickleback has, since the last glaciation, repeatedly colonized and adapted to predator-free freshwater systems, by evolving a low plated phenotype characterized by a reduced number of bony lateral plates (Bell 2001, Leinonen et al. 2011a. This striking example of repeated evolution is characterized by a consistent reuse of globally shared standing genetic variation Conte 2009, Jones et al. 2012). The lateral plate number is a trait under selection (Barrett et al. 2008) and 75% of the variation is determined by a major locus, ectodysplasin (EDA; Colosimo et al. 2004, Le Rouzic et al. 2011. Plates act as armor against gape limited toothed predators (Reimchen 1992, Bergstr€ om 2002, and there is both observational and experimental evidence to support the notion that the low plate phenotype is favored in freshwater populations due to reduced predation by fish (Bell et al. 2004, Cresko et al. 2004, Kitano et al. 2008). Therefore, it is possible that the observed differences in distribution of the different plate phenotypes along the coast in fact reflects a temporal change in predation pressure due to declining predatory fish abundances (see also Jakubavi ci ut_ e et al. 2018).
Although the different stickleback phenotypes co-occur in the Baltic Sea (Jakubavi ci ut _ e et al. 2018), our results indicate that they differ in their diets at least in coastal areas. The diet differences were caused by one prey group only: completely plated sticklebacks more often consumed amphipods. Amphipods are important benthic herbivores in the area (Eriksson et al. 2009, Sieben et al. 2011b, and they form a key grazer group that mediate cascading ecosystem effects in coastal areas across the world (Duffy and Hay 2000, R aberg and Kautsky 2007, Moksnes et al. 2008, Duffy et al. 2015, Donadi et al. 2017. We have no indication as to why completely plated stickleback ate more amphipods than the partially plated individuals. Our expectations were that the partially plated stickleback would seek shelter from predators in the vegetation to a larger degree than the completely plated ones, since they have less armor. We therefore assumed that partially plated stickleback would feed more on those mesograzers that themselves seek shelter and feed in the vegetation, such as amphipods. The opposite may be true if the completely plated stickleback in general are better adapted to habitats with predators; combining better plate armor with being more cautious and seeking shelter to a larger degree. Experimental testing or more controlled sampling of different phenotypes, while controlling for prey availability, is needed both to understand what drives the diet difference and to evaluate how important this diet difference is for food web dynamics. However, we already know that changes in the distribution of stickleback plate phenotypes can have strong effects on ecosystem function in a number of ways (Harmon et al. 2009). Effects of rapid diversification of partially plated stickleback include changes in prey community structure, biomass, and abundance , Rudman and Schluter 2016, Schmid et al. 2019 Fig. 6. Conceptual model illustrating how the density development and trait distribution of a mesopredator population depend on the presence of a higher trophic level of predators and productivity of the system (according to the exploitation ecosystems hypothesis by Oksanen et al. 1981). With top predators, the system has four trophic levels (dotted blue line). Without top predators, the system has three effective trophic levels and mesopredators dominate the predator community (solid red line). (a) In a system with three trophic levels, the mesopredator community accumulates biomass. Ecological effects: Increase in mesopredator density with increasing productivity Oksanen 2000, Eriksson et al. 2012). Evolutionary effects: The absence of toppredators results in a relaxed selection for antipredator traits; and increased selection for traits favorable under density-dependent competition for food (Svanb€ ack and Bolnick 2007). (b) Conceptual model illustrating how increased mesopredator density and resulting changes in trait distribution (caused by the combination of predator loss and increased productivity) may influence trophic dynamics by changing in energy flow (e.g., changes in diet, elemental stoichiometry, behavior), thereby leading to a broad-sense eco-evolutionary feedback on food web structure (Hendry 2017, de Meester et al. 2019). that change the effects of the mesopredator release, experimental studies are needed.
The relative importance of different factors for the increase in three-spine stickleback abundance in the Baltic Sea are currently not known. Our results show that juveniles of the main coastal predator, perch, have decreased over the past decades and that the selection for antipredator traits (plating) may have become relaxed as the incompletely plated phenotype dominates the studied stickleback population. This suggests that decreasing predation from declining abundances of coastal predators has been important for the expansion of the stickleback population in the past decade. Yet, the causality needs to be established in future studies, as there are many other factors that may contribute to the increased stickleback abundance; such as increasing temperatures and habitat change (Eriksson et al. 2011, Lefebure et al. 2014, Nilsson et al. 2014, Sundblad and Bergstrom 2014, Hansen et al. 2019. In fact, similar to the Baltic Sea, the White Sea has shown an increase in three-spined stickleback abundance since the late 1990s (Lajus et al. 2020). The increase in the White Sea has been attributed mainly to increased water temperatures and subsequent changes in habitat availability (Rybkina et al. 2017, Lajus et al. 2020. The migrating behavior of stickleback in the Baltic Sea complicates the analysis of causal factors further by coupling pelagic and coastal processes across scales. The change in the fish community along the coast may therefore be caused by changes both in local and/or pelagic interactions and conditions. Our results, together with evidence for that changes in the trait distribution in three-spined sticklebacks can have variety of community and ecosystem effects (e.g., El-Sabaawi et al. 2016, Rudman and Schluter 2016, Paccard et al. 2018, Schmid et al. 2019), indicate the potential for eco-evolutionary effects in the broad sense (de Meester et al. 2019).
We therefore hypothesize that the changes in the relation between larger predatory fish and their prey, mesopredators, have generate changes in selection regimes, which in turn may feedback on a variety of ecological processes (Fig. 6). Models of trophic dynamics suggest that predator effects on lower trophic levels interact with resource availability, such that the abundance of the highest (apex) trophic level increases with primary production (Fretwell 1977, Oksanen et al. 1981, Nisbet et al. 1997, Oksanen and Oksanen 2000, Eriksson et al. 2012. Thus, in a system where larger predatory fish control trophic dynamics, the larger predator community accumulates biomass. For their mesopredator prey, there should be an increase in selection for antipredator traits, but no increase in density with increasing productivity Oksanen 2000, Eriksson et al. 2012). However, in a system where the large predators have declined strongly, energy (biomass) accumulates at the mesopredator level, while at the same time the absence of top predators result in a relaxed selection for antipredator traits. Here, increasing productivity would lead to an increase in mesopredator density, and consequently, densitydependent competition (Oksanen and Oksanen 2000, Svanb€ ack and Bolnick 2007, Eriksson et al. 2012. Thus, the abundance of larger predators and the productivity of the ecosystem should together determine mesopredator abundance as well as trait distribution (Fig. 6a), which may in turn influence trophic dynamics through changes in, for example, diet composition, behavior, and/ or elemental stoichiometry; and thereby generate an eco-evolutionary feedback on food web structure (broad-sense eco-evolutionary feedback; Hendry 2017, de Meester et al. 2019, Fig. 6b).
The hypothesis of an eco-evolutionary feedback in the Baltic Sea is partly supported by our results: The number and proportion of incompletely plated stickleback increased with production in bays with low abundances of perch, indicating a reduced selection for antipredator traits and increased density-dependent selection when larger predators of stickleback decrease. We also found small but significant differences in diets between the phenotypes, where the incompletely plated individuals consumed less amphipods; which is a key benthic herbivore regulating trophic transfer between predators and algae in many shallow marine systems across the world (e.g., Duffy and Hay 2000, Moksnes et al. 2008, Duffy et al. 2015, Donadi et al. 2017. However, to provide solid empirical predictions of eco-evolutionary consequences of the shift in selection pressures, we need experimental testing of the ecological effects of different phenotypes, behavioral studies to assess v www.esajournals.org morphology-dependent mating patterns, and genetic studies to identify signatures of selection. Our current results emphasize the need to better understand how ecological feedbacks can change selective regimes, potentially driving rapid adaptation, and the need to take these effects into account in ecosystem management strategies (de Meester et al. 2019, Hendry 2019).

ACKNOWLEDGMENTS
We thank the property and fishing right owners around each of the 32 study bays for facilitating the fieldwork and a number of students, volunteers, and fish biologists for assistance in the field and laboratory: E. Anderberg, F. Ek, G. Johansson, P. Jacobson, C. J€ onander, E. M€ ork, L. N€ aslund, O. Pettersson, G. Lil-liesk€ old Sj€ o€ o, S. Skoglund, M. van Regteren, M. van der Snoek, V. Thunell and L. Vik. We thank Veijo Jormalainen and Joel Fodrie for helpful comments on an earlier version of the manuscript. The authors declare no conflict of interest. U. Bergstr€ om and J. S. Ekl€ of contributed equally to the work reported here.