Two closely related species differ in their regional genetic differentiation despite admixing

Genetic differences between regions are usually studied for individual species. However, many species can reproduce with each other. We studied whether gene flow between two closely related sedge species influences regional differences. Our molecular genetic data support considerable gene flow between the species. Still, we detected clear genetic differences between species and regions, and more pronounced regional differences for the less common one. Thus, gene flow between the species appeared too weak to neutralize differences between the regional genetic structure of our study species. We encourage further regional differentiation studies in groups of cross-compatible species.


Introduction
Plants and other organisms differ in their levels of genetic diversity and genetic differentiation (Linhart and Grant 1996). The extent of genetic differentiation among populations and regions depends on the balance of evolutionary forces decreasing and increasing genetic differentiation, that is, gene flow, genetic drift, mutation and selection (Slatkin 1987). The relative importance of these forces may be affected by selection strength, as well as population size, environmental barriers to dispersal and plant life history traits, especially mating system and dispersal mechanism (Loveless and Hamrick 1984). Higher differentiation among populations is generally found for clonally reproducing and selfing species (Loveless and Hamrick 1984;Hamrick and Godt 1996;Nybom 2004;Song et al. 2006) and for species with disjunct distributions or small populations (Ellstrand and Elam 1993), due to effects of genetic drift or reduced gene flow (Slatkin 1987;Schönswetter et al. 2006). Outcrossing and sexually reproducing species, conversely, show less differentiation among populations (Hamrick et al. 1979;Loveless and Hamrick 1984).
While genetic consequences of small population size, i.e. increased inbreeding and genetic drift, are expected to contribute to higher population differentiation and lower genetic diversity within populations (Slatkin 1987;Ellstrand and Elam 1993;Amos and Harwood 1998;Premoli 2003;Leimu et al. 2006;Zhivotovsky et al. 2016), this is not found in all studies (Gitzendanner and Soltis 2000). Loveless and Hamrick (1984) agree that small populations are more susceptible to drift and fixation, but suggest that immigration is more effective in altering gene frequencies in small populations. Thus, gene flow can prevent differentiation and loss of genetic variability, especially in long-lived species. Levels of genetic diversity in populations may also depend on distance from glacial refugia (Schönswetter et al. 2006) and on the landscape (Holderegger and Wagner 2006). For populations at higher altitudes higher radiation intensity may increase the rate of mutations and thus population genetic diversity (Li et al. 1997).
It is generally assumed that species are reproductively isolated without gene flow between them, but in reality hybridization is a widespread phenomenon (Ellstrand et al. 1996;Abbott et al. 2013) and becomes more frequent with continuous climate change and human influence on the environment (Hoffmann and Sgro 2011). Gene flow between taxa can have a profound effect on genetic diversity (Derieg et al. 2008), and due to subsequent genetic drift, may lead to population differentiation (Bain and Golden 2003). Hybridization of closely related species is found to lead to genetic differentiation among regions (Rieseberg et al. 1999;Kane et al. 2009;Lepais et al. 2009;Krebs et al. 2010;Brennan et al. 2016). According to the 'semipermeable species boundaries' theory, alleles at some loci can be exchanged between species and species boundaries can vary geographically (Harrison and Larson 2014). Thus, hybridization coupled with backcrossing can be expected to affect population differentiation on a regional scale.
Hybrid individuals are expected to be more heterozygous than their parental taxa due to genetic admixing (Rieseberg and Wendel 1993;Allendorf et al. 2001;Barton 2001;Harrison and Larson 2014;Todesco et al. 2016). Hybridization could result in beneficial evolution by providing additional adaptive genetic variation (Lewontin and Birch 1966;Arnold 2006;Abbott et al. 2013). Conversely, hybridization can reduce biodiversity by causing loss of alleles and genetic diversity by genetic assimilation (Levin et al. 1996). Genetic studies on hybridizing taxa have focused on their genetic variation or the extent of hybridization and introgression (e.g. Friedman et al. 2008;Volkova et al. 2008;Kane et al. 2009;Korpelainen et al. 2010;Krebs et al. 2010). Our study addresses the differentiation of hybrids from parental taxa.
We investigate genetic diversity and regional differentiation between and within Carex flava, C. viridula and their hybrid C. × subviridula for populations from three regions (Estonia in Northern Europe and Lowland and Highland Switzerland in Central Europe). Carex flava and C. viridula var. viridula sensu stricto (s.s.) (henceforth C. viridula) are wind-pollinated, self-compatible, caespitose perennials of the C. flava aggregate (Carex sect. Ceratocystis, Cyperaceae). Although there are no impediments to outcrossing, a large amount of seeds is produced by selfing (Vonk 1979;Schmid 1984a). Both taxa, with circumpolar distribution, occur in the temperate and subarctic Northern hemisphere and also in North Africa (Hultén and Fries 1986;Crins and Ball 1989;Koopman 2011). They often co-occur and hybridize, especially at sites with disturbances (Vonk 1979;Jiménez-Mejías et al. 2012), resulting in C. × subviridula. Carex viridula is considered to be a dispersal generalist, potentially being transported by biotic, e.g. birds, mammals, invertebrates, and abiotic agents, e.g. water and wind (Schmid 1984a;Crins and Ball 1989). Long-distance dispersal is proven experimentally with C. flava var. alpina seeds remaining still partly intact after passing the digestive tract of domesticated ducks within 18 h (Schmid 1984a). Carex viridula is a weak competitor, but able to colonize and survive in fluctuating, relatively unpredictable moist or wet habitats, where it forms small populations (Schmid 1986;Crins and Ball 1989;Kuchel and Bruederle 2000). However, populations of C. viridula are sensitive to anthropogenic influence, such as the drainage of mires, regulation of water levels and eutrophication of shores (Pykälä and Toivonen 1994) and its decrease in southern parts of its distribution could be explained by climate change (Parmesan and Yohe 2003;Chen et al. 2011). Hence, the occurrence of C. viridula has decreased over time and populations have become more fragmented (Davies 1953;Pykälä and Toivonen 1994). According to the latest red list of endangered plants in Switzerland, it is considered near threatened due to loss of habitat (Bornand et al. 2016). Carex flava, on the other hand, is a strong competitor and not as sensitive to environmental changes; its populations are larger and rather constant in time (Schmid 1984a(Schmid , b, 1986. Intraspecific gene flow is expected to be higher for C. flava and smaller for C. viridula (Schmid 1986).
Previous genetic studies with the C. flava agg. were restricted to small regions, did not consider co-occurrence and hybridization or used allozyme markers of limited variability (Bruederle and Jensen 1991;Kuchel and Bruederle 2000;Hedrén 2002Hedrén , 2004Blackstock and Ashton 2010;except Jiménez-Mejías et al. 2012). The novelty of our study originates from considering co-occurrence and hybridization, from comparing geographically and climatically distant regions, and from using contemporary molecular microsatellite (SSR) markers. According to our hypothesis, admixture, i.e. hybridization coupled with backcrossing, can be expected to affect population differentiation on a regional scale. Therefore, species characterized by relatively small populations might be less well differentiated on a regional level than expected without admixing. Thus, we study whether, despite admixing, C. viridula, characterized by relatively small populations, exhibits higher levels of population differentiation than C. flava, which has relatively large populations. In addition, we investigate the variability of hybrid populations and their differentiation from the parental taxa.

Study taxa and regions
In the studied regions, the C. flava group comprises the four taxa C. flava, C. lepidocarpa, C. demissa and C. viridula var. viridula (Schmid 1981;Toom et al. 2016). In Estonia, the two varieties C. viridula var. pulchella and var. bergrothii are also found (Toom et al. 2016). In this study, we follow Hedrén (2002Hedrén ( , 2004) taxonomic treatment, but we use the more common name C. viridula instead of calling it C. oederi.
The studied regions differ in their postglacial history. Swiss populations were established earlier after glaciation, because the territory of Estonia was covered by ice at the end of the last ice age 18 000 BP, when the lower parts of Switzerland were ice-free (Hewitt 1999). In addition, the current populations in Switzerland are much closer to possible southern glacial refugia in Iberia, the Apennine or the Balkan Peninsula (Schönswetter et al. 2006) than Estonian populations are. As the populations in Switzerland are older, they may be more amalgamated via interspecific gene flow, hybridization and introgression than the Estonian ones.

Population sampling
We collected 380 samples of C. flava, C. viridula and C. × subviridula populations from seven sites in Highland Switzerland (>1000 m, in 2012), five sites in Lowland Switzerland (<1000 m, in 2012) and 12 sites in Estonia (2013; Table 1; Fig. 1). At 15 sites we found C. flava and C. viridula growing together or with other sedges, with whom they are able to hybridize, i.e. with other C. flava aggregate members or C. punctata and C. hostiana (Davies 1955;Schmid 1982;Crins and Ball 1989;Więcław and Wilhelm 2014). Samples of the C. flava × C. viridula hybrid (C. × subviridula) were found at two sites in Estonia and two in Lowland Switzerland, while none were found in Highland Switzerland. As hybrids, we classified partly sterile individuals that were morphologically intermediate between the parents, but often were more robust and pale (Schmidt et al. 2017). The sampling populations of C. viridula were smaller, i.e. had fewer individuals, than the ones of C. flava. As we focus on between-region and between-taxa comparisons, for which the population is the unit of replication against which between-region and between-taxa differences are tested, and as the power of such analyses depends on the number of populations, whereas the replication within populations is less decisive (van Kleunen et al. 2014), we sampled a total of 37 populations and on average 10.3 individuals per population.
As population genetic diversity is expected to be high in the centre and to decline at the margins of distribution ranges (Volis et al. 2016), we studied the taxa neither in the centre nor the margins of theirs. Our southernmost population was at the latitude of 45.96 N in Caslano, Switzerland (V36), and the northernmost one at 59.29 N in Anija, Estonia (CF13). In Switzerland, the population at lowest altitude was in Caslano at 270 m a.s.l. and the one at highest altitude in Arosa at 1920 m a.s.l. Vouchers with samples of all study populations were deposited in the herbarium of the Natural History Museum of the University of Tartu.

Microsatellite analysis
The microsatellite loci tested in this study were developed for other species. Due to good success in crossamplification among congeners (Rossetto et al. 1999), we chose primer pairs isolated from other Carex species. We chose in total 17 polymorphic microsatellite loci from previous studies that had showed successful cross-species amplification. We screened nine primer pairs developed for Carex scoparia (Hipp et al. 2009), two primers developed for Carex rugulosa (Ohbayashi et al. 2008), four primers developed for Carex kobomugi (Ohsako and Yamane 2007) and two primers developed for Carex limosa (Escudero et al. 2010). Primary microsatellite analysis was performed with few individuals of each taxon and 17 primer pairs. Of the 17 primer pairs tested, 10 aligned successfully with recipient DNA, cross-amplified in C. flava, C. viridula and C. × subviridula, exhibited polymorphism and showed identifiable peaks in fragment analysis. Those 10 primer pairs were used for further analysis with all 380 samples [see Supporting Information- Table S1]. We did not sequence the fragments recovered for the 10 loci nor did we perform progeny analysis. However, the primer pairs used in our study had been successfully used in other population genetic studies (Hipp et al. 2009;Escudero et al. 2010;Korpelainen et al. 2010). Each primer was optimized for a range of temperatures (T a : 49-60.1 °C). Magnesium source 1.2 mM MgSO 4 was used, except for S177, where 1.6 mM MgCl 2 was used. PCR amplifications were performed in 10 µL volumes containing 1-2 µL of genomic DNA, 1.2 µL GoTaq Flexi buffer (1×), 0.6 µL of each dNTP, 0.5 µL of untagged primer, 0.5 µL of fluorescent tag, 0.5 µL of the tagged primer, 0.05 µL of bovine serum albumin (BSA), 0.05 µL GoTaq Flexi DNA polymerase and varying concentrations of MgCl 2 or MgSO 4 . Total genomic DNA was isolated from silica-dried leaves using the CTAB method (Doyle 1987). The extracted DNA was dissolved in 100 µL of TE buffer and diluted to 1:10 for further PCR analyses. DNA from each sample was amplified with the common tag containing one of four fluorescent dyes, 6-FAM, PET, VIC or NED (Applied Biosystems). PCRs were carried out as follows: preliminary denaturation at 95 °C for 5 min, 35 cycles at 95 °C for 1 min, annealing temperature 53.6-60 °C for 1 min, 72 °C for 1 min and a final extension step at 72 °C for 30 min, using a Techne TC-5000 thermocycler (Bibby Scientific). PCR products of different primers, each of 1-2 µL, were mixed together yielding a total of 20-µL mixture. From each mixture 2 µL were pooled with 10 µL buffer (size standard: deionized formamide = 1:25) in the wells of a 96-well plate for fragment analysis on a 3730xl DNA Analyzer (Applied Biosystems), where samples of different plant individuals were randomized. Product sizes were determined using the Peak Scanner Software v1.0 (Applied Biosystems). Scoring errors, e.g. null alleles, were identified and corrected using micro-checker 2.2.3 (van Oosterhout et al. 2004).

Data analysis
For each population the effective number of alleles (N e ), percentage of polymorphic loci (PL), expected heterozygosity (H e , also called gene diversity) and observed heterozygosity (H o ) were estimated across all loci using GENALEX 6.5 (Peakall and Smouse 2006). Allelic richness (Ar) and compliance with Hardy-Weinberg expectations were calculated in FSTAT v 2.9.3 (Goudet 2002). The inbreeding coefficient (F IS ) shows the probability to observe alleles of an individual (I) that are identical by descent (IBD) in a subpopulation (S). F IS , calculated as (H e − H o )/H e , allowed to estimate the prevailing mating systems by region (i.e. Estonia, Lowland and Highland Switzerland) and  Table 1. taxon. Variation in these measures of population diversity was tested with ANOVA, using taxon, region of origin and the interaction of taxon and region of origin as fixed effects, implemented in the software R v 3.1.2 (R Development Core Team 2016).

Population genetic structure and hybrid identification.
We used a Bayesian clustering approach as implemented in STRUCTURE v 2.3.4 (Pritchard et al. 2000) (i) to estimate the number of genetic clusters (K) without a priori knowledge of taxonomy or population, and (ii) to identify the hybrid individuals with admixture analysis. The clustering was conducted with the admixture model and the correlated-allele-frequencies option using a burn-in of 10 000 steps and 100 000 replications, the remaining parameters were set to the default values. Five independent runs were done for the set of K = {1:10}. We used Structure Harvester v 0.6.92 (Earl 2012) to visualize the optimal number of clusters (K) by using firstly the ΔK method of Evanno et al. (2005) and secondly by examining the distribution of the log-likelihoods for the value with the highest probability and lowest variance using Markov chain Monte Carlo simulations. With K = 2 we detected the posterior probability (q-values), which describes the proportion of an individual genotype originating from each of K categories. We used K = 2 as we expected two taxa contributing to the gene pool of hybrids. We chose a threshold value of 0.9, which was found efficient to distinguish pure individuals (q > 0.9 or q < 0.1) from hybrids and backcrosses (0.1 < q < 0.9) (Vähä and Primmer 2006;Burgarella et al. 2009).

Genetic differentiation.
To compare the degree of differentiation among groups of populations categorized by taxa and region of origin, between-group F ST values were calculated, using Arlequin v 3.5.2.2 (Excoffier et al. 2005). The significance of differences in the resulting values was tested with 1000 permutations. To illustrate the dissimilarities among groups of populations categorized by taxa and by region of origin, multidimensional scaling (MDS) analysis was performed in software R v 3.1.2 (R Development Core Team 2016) based on Reynolds distances obtained with the software Arlequin, which estimates the co-ancestry of different samples (Reynolds et al. 1983).
Analysis of molecular variance (AMOVA) enabled us to determine the distribution of microsatellite variation among groups of populations, among populations within groups and among individuals within populations, using Arlequin. We grouped the populations according to the tested hypotheses per taxa and regions. We tested the significance of the variance components by calculating their probabilities based on 9999 permutations of individual samples.

Results
The total number of alleles observed per locus in the overall sample of 380 individuals from 37 populations and 24 sites of three regions ranged from 4 to 11, with overall 64 alleles scored over the 10 loci. Private alleles were detected at seven loci, four of which were found in C. flava in Highland Switzerland [see Supporting Information- Table S2].
Total genetic diversity varied little between the three studied taxa, and the overall absolute values of genetic diversity statistics, which comprise variation within and between regions, were similar in the more common C. flava (N e = 1.42, PL = 56.11 %, Ar = 1.56, H e = 0.21) than in the less common C. viridula (N e = 1.54, PL = 54.0 %, Ar = 1.64, H e = 0.25) ( Table 2). The percentage of polymorphic loci was significantly different between two regions (F = 10.2, P = 0.01) and among three regions (F = 5.85, P = 0.01; Table 3). Significant taxon-by-region interactions for allelic richness (F = 3.58, P = 0.04), percentage of polymorphic loci (F = 5.30, P = 0.01) and expected heterozygosity (F = 3.89, P = 0.03) indicate that differences between the two taxa in their levels of genetic diversity depended on the region (Table 3a).

Population diversity and inbreeding
The genetic diversity in C. viridula populations was highest in Lowland Switzerland (N e = 1.77, PL = 70.0 %, Ar = 1.86, H e = 0.32), followed by Estonia (N e = 1.49, PL = 51.1 %, Ar = 1.59, H e = 0.23). In Highland Switzerland only few populations of C. viridula were found and their genetic diversity was much lower than in the other two study regions (N e = 1.33, PL = 35.0 %, Ar = 1.40, H e = 0.16). We expected C. viridula with its smaller populations to be more prone to inbreeding. Accordingly, C. viridula showed a deficit of heterozygotes in Lowland Switzerland populations and in some populations of Estonia. However, other Estonian populations and the    (Table 2b) (Table 2a), the heterozygote excess was not significant, however. These results could be explained first with proximity to several glacial refugia and putative mutagenic influence of higher radiation in mountains.
Genetic diversity was similar for hybrids in Lowland Switzerland (N e = 1.49, PL = 65.0 %, Ar = 1.68, H e = 0.25) and in Estonia (N e = 1.58, PL = 45.0 %, Ar = 1.69, H e = 0.24). In both regions hybrid populations showed an excess of heterozygotes (Table 2c), and were more variable than either parental taxon, in line with expectations for populations with genetic admixing.

Nuclear admixture analysis
The ΔK method in the Bayesian program STRUCTURE classified all individuals into two clusters [see Supporting Information- Fig. S1], illustrating a clear distinction between taxa, but not among regions within the taxa (K = 2; Fig. 2A). We used the admixture coefficients from the analysis with K = 2 to determine the proportions of admixed individuals. The majority of C. flava and C. viridula individuals from 'pure' populations were indeed classified with high admixture coefficients (q > 0.90), i.e. as pure individuals, only two putative C. flava and two C. viridula individuals had q < 0.90, suggesting that they rather were mixed genotypes.
The situation was more complex at 'mixed' sites, where two or more taxa of the C. flava complex cooccurred. For the C. flava morphotypes at the mixed sites 136 individuals of 159 putative C. flava were indeed pure C. flava (85.5 %; q > 0.90), 21 of 159 had mixed genotypes (13.3 %; q < 0.90) and two individuals (1.2 %) were even C. viridula-like. For the C. viridula morphotypes at the mixed sites, 102 individuals of 110 putative C. viridula were pure C. viridula (92.7 %; q > 0.90), three had mixed genotype (2.7 %; q < 0.90) and five (4.5 %) were C. flava-like. For the C. × subviridula hybrid morphotypes a wide range of admixture proportions were found (q ranged from 0.103 to 0.897), suggesting the presence of a broad range of hybrid generations and backcrossing to both parental taxa.
With the log-likelihood distribution method, the value where the rate of increase in likelihood reaches a plateau without increase in variance corresponds to three clusters [see Supporting Information- Fig. S1]. K = 3 yielded different results than K = 2, also revealing differences between regions, namely admixture between plants of C. flava and C. viridula in Highland and Lowland Switzerland (Fig. 2B). The third cluster (blue cluster with q > 0.1 in Fig. 2B) occurred mostly for individuals at  Table 1. mixed sites of C. flava and C. viridula (in seven of nine of the mixed sites, and also in Kaltbrunn (F33) and in Melchsee-Frutt (F1), where we exclusively found C. flava).
Overall, the STRUCTURE results indicate clear differentiation between the taxa and in addition further variation between regions. This is in line with the results of genetic diversity, which showed significant differences between regions (for the percentage of polymorphic loci) and significant region-by-taxon interactions (for three measures of genetic diversity; Table 3).

Interspecific differentiation
According to the hierarchical AMOVA, the proportion of genetic variance within regions between the studied three taxa was highest in Estonia (39.16 %), followed by Highland (28.27 %) and Lowland Switzerland (17.51 %), i.e. between-taxa differentiation was highest in Estonia and lowest in Lowland Switzerland, where three of the five studied sites were mixed (Table 4)
Ordination according to the MDS analysis illustrated the genetic distances of populations grouped by taxonomic identity and region of origin (Fig. 4). Estonian C. flava populations formed a distinct cluster, whose difference from clusters of Lowland and Highland Switzerland C. flava populations was smaller than the differences observed among C. viridula populations between the regions (Fig. 4). Estonian C. viridula populations clearly differed from Lowland Switzerland populations and even more from the two populations of C. viridula in Highland Switzerland (V9 and V12), which were also very different from each other (Fig. 4). These findings clearly suggest higher differentiation in C. viridula than in C. flava.
Moreover, in combination, the findings on interspecific and intraspecific differentiation indicate that differentiation between the taxa was stronger than differentiation within the taxa between regions. This was further supported by hierarchical AMOVA, where 26.56 % of the variation resided between the three taxa, while 13.00 % resided between regions (Table 4a and b). As expected, the AMOVA analyses without hybrids C. × subviridula showed higher variance between the parental taxa (32.02 %; Table 4c).

Discussion
Our microsatellite results support growing evidence that interspecific gene flow is more widespread than previously suspected, but that between-species differences are still retained by various mechanisms (Arnold et al. 1990;Friedman et al. 2008;Smith and Waterway 2008;Kane et al. 2009;Nolte et al. 2009;Scascitelli et al. 2010). Microsatellite markers are proven appropriate for population structure and differentiation studies (Pálsson 2000;Song et al. 2006;Tyagi et al. 2016;Volis et al. 2016) and for investigating the relationship among closely related taxa (Korpelainen et al. 2010;Talve et al. 2013Talve et al. , 2014. We examined whether C. viridula, the taxon with a more disjunct distribution and smaller populations, showed lower genetic diversity and higher inbreeding than the more widespread C. flava, whose populations are larger and more constant in time. However, mean gene diversity (expected heterozygosity) was slightly higher in C. viridula, though not significantly higher from C. flava with H e = 0.25 and H e = 0.21, respectively. Carex viridula showed highest genetic diversity in Lowland Switzerland, whereas C. flava was most diverse in Highland Switzerland, with H e = 0.32 and H e = 0.26, respectively. Both taxa had lower diversity in Estonia. Greater allozyme diversity and a lower inbreeding coefficient for C. viridula than for C. flava was also detected in earlier studies using allozymes by Hedrén (2004) and Bruederle and Jensen (1991), where the latter had considered C. viridula s.l., united with C. lepidocarpa, C. demissa and C. viridula s.s., however. Meanwhile, Kuchel and Bruederle (2000) detected low levels of allozyme diversity in North American C. viridula and attributed it to bottlenecks at arrival from Europe and to predominant selfing. Higher diversity in highland populations, as in our study in C. flava, was found in some studies, e.g. for Cystopteris fragilis (Pteridophyta) using isozymes (Gämperle and Schneller 2002), for Primula farinosa (Primulaceae) using RAPD analysis (Reisch et al. 2005) and for Campanula thyrsoides (Campanulaceae) (Frei et al. 2012), and was suggested to be due to higher mutation rates due to elevated radiation (Li et al. 1997). However, others reported higher diversity of lowland populations, as we found for C. viridula (e.g. Premoli 2003;Schönswetter et al. 2006).
Deviation from HWE may indicate inbreeding. In earlier studies, high selfing and evidence for inbreeding has been found in Carex (Arens et al. 2005;King and Roalson 2009;Escudero et al. 2010;Kull and Oja 2010). Inbreeding has also been found to predominate in selfcompatible caespitose sedges . We found an excess of heterozygotes in C. flava when Table 4. AMOVA of Carex flava (sect. Ceratocystis, Cyperaceae), C. viridula var. viridula and C. × subviridula for SSR data considering the whole data set of all 37 populations with two or three hierarchical levels (a and b) or subsets of populations with two or three hierarchical levels (c-i). P-value = associated significance derived from 16 000 permutations. *P < 0.001; **P < 0.05.

Data set
Source of variation df % of variation  (Ohsako and Yamane 2007;Hipp et al. 2009). Earlier it has been shown that C. flava is the main partner for backcrossing, as it can occasionally be pollinated successfully by F 1 hybrids or backcrosses (Schmid 1982;Więcław and Wilhelm 2014). This suggests that C. flava was more prone to between-taxa crosses when growing adjacent to other taxa in section Ceratocystis, e.g. C. viridula, C. lepidocarpa and C. demissa. On the other hand, the direction of hybridization may be affected by the length of the style (Field et al. 2011) which in carices would be determined by the beak length, and hybridization may occur more commonly from long-beaked carices to short-beaked carices due to pollen competition (Derieg et al. 2013), which would suggest higher gene flow from C. flava to C. viridula than vice versa. We observed hybrids in most of the sites where C. flava and C. viridula grew together sympatrically. Mean gene diversity (expected heterozygosity) was not higher in hybrid populations than in the parental taxa (H e = 0.24 vs. 0.21 and 0.25). Korpelainen et al. (2010) studied sedge hybrids in Finland using microsatellite data and found for C. aquatilis × recta (sect. Phacocystis) similar gene diversity than for its parental taxa (H e = 0.348 vs. 0.308 and 0.460) and for C. paleacea × recta higher diversity (H e = 0.603 vs. 0.185 and 0.460). In contrast, high genetic diversity, using RAPD analysis, was found in Fallopia × bohemica (Polygonaceae) in Germany and Switzerland (Krebs et al. 2010).
We determined hybrids based on partial sterility and morphological differences from the parental taxa. Most hybrid individuals had utricles without fully developed achenes, but some hybrid individuals had circa 5 % of fully developed achenes (own observation). The admixture proportions detected by microsatellites for hybrid individuals were very variable, indicating that these comprised F 1 to F n hybrids and backcrosses. Our genetic data showed that some of the supposed intermediate individuals were not hybrids, but rather backcrosses, or in rare cases even pure parental taxa. Thus, our results imply that hybrids can be identified well based on morphological criteria, but that morphological criteria do not allow for distinguishing backcrosses and that they may   Table 1. even lead to occasional, but very rare, errors in identifying pure individuals.

Interspecific differentiation
The results of our STRUCTURE analysis suggested that we dealt with two species and their hybrids (K = 2) and with some regional differentiation (K = 3; Fig. 2). The presence of separate taxa was further supported by hierarchical AMOVA, where we detected 13.00 % of genetic variation among studied regions, but more than twice this variation between taxa (32.02 %; Table 4). We conclude that there are solid differences between the two species despite evident hybridization. This is in line with earlier studies on the taxonomic relationships of C. flava and C. viridula (Bruederle and Jensen 1991;Hedrén 2002;Jiménez-Mejías et al. 2012). Similarly, Morgan-Richards and Wolff (1999) found differences between sympatric Plantago major (Plantaginaceae) taxa preserved despite intraspecific gene flow. Kane et al. (2009) also concluded that hybridizing Helianthus taxa (Asteraceae) remained largely reproductively isolated and morphologically and ecologically distinct despite high levels of interspecific gene flow. Brennan et al. (2016) found similar results for hybridizing Senecio taxa (Asteraceae) and explained them with selection against hybrids and locally maladapted hybrid individuals.
Our STRUCTURE analysis with K = 3 revealed a widely present third genetic cluster in Switzerland, which is extremely rare in Estonia (Fig. 2). This third cluster occurred in all three studied taxa C. flava, C. viridula and C. × subviridula, at both low and high altitudes in Switzerland. This suggests high gene flow between taxa within Switzerland and lower between Estonia and Switzerland, as further supported by AMOVA (Table 4). In Estonia, hybridization and introgression occurs, whereas in Lowland Switzerland gene flow between the species seems to be more frequent and to affect the genetic structure of C. flava and C. viridula.

Intraspecific differentiation between regions
As expected for rarer taxa with smaller populations, we found higher among-region differentiation for C. viridula (AMOVA; 20.77 %) than for the more widespread C. flava or hybrid (6.84 and 18.27 %, respectively). In addition, groupwise F ST values between regions were higher for C. viridula, especially between Estonia and both altitudes in Switzerland (Fig. 3), indicating that hybridization was not strong enough to prevent stronger regional differentiation in C. viridula. High levels of differentiation might be caused by low levels of wind pollination, which is not very effective for small herbs in closed habitats (Kull and Oja 2010). Another explanation for the higher differentiation of C. viridula is the loss of suitable habitats and subsequent fragmentation of populations, as also shown for other taxa (Pykälä and Toivonen 1994;Reisch et al. 2005;Schönswetter et al. 2006;Kull and Oja 2010;DeWoody et al. 2015). Taking into account that C. viridula had a wider distribution in the past and now shows reduced occurrence and increased fragmentation, our finding of higher differentiation between regions for C. viridula fits theoretical expectations.
We found especially low genetic differentiation among the C. flava groups of populations between Highland and Lowland Switzerland (Fig. 3). This supports the idea of higher gene flow between populations of C. flava, as suggested by Schmid (1984bSchmid ( , 1986. With the use of microsatellite markers, Escudero et al. (2013) have shown considerable mixing among populations in the widespread C. scoparia (sect. Ovales) and explained it with long-distance dispersal. Potential for long-distance pollen and seed dispersal was suggested to contribute to low geographic differentiation of circumpolar C. bigelowii (sect. Phacocystis) (Schönswetter et al. 2008). Theoretical studies have shown that only a small amount of long-distance gene flow is needed to prevent population differentiation for neutral alleles (Loveless and Hamrick 1984). Differences in phenology with altitude are expected to reduce gene flow via pollen and increase differentiation instead (Premoli 2003;Reisch et al. 2005). As flowering times differ notably in Switzerland (Körner 1999), we suggest that this may explain the observed low differentiation among Swiss C. flava populations of similar altitude.
Hybrids showed unexpectedly high differentiation between Estonia and Switzerland (F ST = 0.32). This difference could originate from genetic drift in the hybrid populations (Nolte et al. 2009). Moreover, backcross patterns may have differed between the regions. A further explanation could be that in natural populations hybrid swarms involve more than two species (Lepais et al. 2009). At some sites, we found C. lepidocarpa and C. demissa growing beside C. flava and C. viridula, which might increase the number of potential parental taxa and differentiation of hybrids. These mechanisms of hybrid and backcross differentiation between regions are very different from the case of hybrids in the invasive Fallopia species complex, which arise from crosses between different taxa in their home origin, and where different hybrids were introduced to different regions, leading to high regional differentiation in the introduced range (Krebs et al. 2010).

Conclusions
Our in-depth analysis of 380 individuals belonging to two sedge taxa and their hybrids and involving populations from three regions suggest that hybridization and introgression are neither strong enough to prevent clear differentiation between taxa nor to prevent stronger regional differentiation for the less common taxon. We encourage further studies on regional differentiation of hybrids and parental taxa to see whether our findings for the C. flava complex represent a more general pattern. Moreover, we suggest also considering hybrids and closely related taxa when addressing genetic diversity and differentiation for rare and endangered taxa.

Contributions by the Authors
L.S., T.O. and M.F. conceived the study and its design. L.S. conducted field and lab work, analysed data and led the writing of the manuscript. T.O. advised lab work. Both T.O. and M.F. contributed to data analysis, data interpretation and writing of the manuscript. This manuscript forms part of the PhD thesis of the first author (L.S.). All authors have made a substantial contribution to the paper.