Genetic Diversity of Stratiotes aloides L. (Hydrocharitaceae) Stands across Europe

Intense land use and river regulations have led to the destruction of wetland habitats in the past 150 years. One plant that is affected by the reduction in appropriate habitats is the macrophyte Stratiotes aloides which has become rare in several areas. The preservation of genetic diversity within a species is a prerequisite for survival under changing environmental conditions. To evaluate the level of genetic diversity within and among populations of Stratiotes aloides, we investigated samples from waterbodies across Europe using AFLP. Low genetic diversity among samples from the same population was found, proving that stands consist of few clones which propagate clonally. Nevertheless, most populations showed differences compared to other populations indicating that there is genetic diversity within the species. The analyzed samples formed two groups in STRUCTURE analyses. The two groups can be further subdivided and mainly follow the major river systems. For conserving the genetic diversity of Stratiotes aloides, it would thus be preferable to focus on conserving individuals from many different populations rather than conserving selected populations with a higher number of individuals per population. For reintroductions, samples from the same river system could serve as founder individuals.


Introduction
The monotypic genus Stratiotes includes the sole living species S. aloides L. (water soldier) and is a member of the Hydrocharitaceae which belong to the order Alismatales within monocots. During the Tertiary and Quaternary periods, there were up to twenty different species of the genus Stratiotes in Europe and Asia [1] (and references therein). The free-floating aquatic macrophyte is perennial with leaves up to 40 cm long and 4 cm wide which are arranged in rosettes. Depending on the season, the plants are emerged or submerged [2]. During the vegetative and reproductive period of a year, the plants are mostly emergent with the rhizoids free in the water or loosely attached to the soil. In autumn, the plant submerges in order to overwinter at the bottom of the water until the following spring [2]. Besides sexual reproduction, the dioecious plants also propagate via vegetative organs (turions and offshoots). Since vegetative reproduction is much more common in S. aloides, stands in one waterbody are often formed by clonal individuals of the same sex [2,3]. As long as individuals from different sexes are not transferred from one waterbody to another by floods and high waters, sexual reproduction is very rare. As a typical flood plain species, it inhabits slow-moving or stagnant waters such as ponds, canals, ditches and oxbow waters where it often dominates macrophyte communities [1]. Stands of water soldiers are frequently inhabited by macroarthropod fauna of which several species are of conservation concern [4][5][6]. Stratiotes aloides is distributed from northern Middle Europe in the West to Siberia in the East.
Results based on uncorrected p-distances and Hamming distances gave the same clustering patterns in neighbor-joining (NJ) dendrograms and principal coordinate analysis (PCoA). The same was found for the pair of Dice distances and Jaccard distances. Therefore, we used only results based on uncorrected p-distances and Dice distances for further analyses.

Neighbor Joining
NJ dendrograms based on uncorrected p and Dice distances both showed a star-like shape with a backbone of relative short branches lacking bootstrap support greater than 75% ( Figure 1). They differed slightly in clustering patterns, but all of the differing branching patterns did not receive high bootstrap support in either of the two analyses. The groups found in STRUCTURE analyses and in PCoA are partly found in the NJ dendrograms. The two groups "Baltic + Hungary" (BH) and "Central European Highlands and plains 1 + Romania" (CER) based on STRUCTURE analyses (K = 5) are supported with high bootstrap values in the NJ dendrograms (BH: 99.7% in Dice, 89.8% in uncorrected p; CER: 99.9% in Dice and uncorrected p). Here, we present only unrooted trees due to the low resolution of their backbone.

STRUCTURE
STRUCTURE analysis gave the highest value of ΔK for K = 3 plus a few other suboptimal K values ( Figure S1a) in the analysis of the reduced dataset (one or two representative individuals per population). However, the latter contained clusters with negligible membership ("empty" clusters). Visualization of K = 45 based on the STRUCTURE results (reduced dataset) showed six clusters which are subsets of the

STRUCTURE
STRUCTURE analysis gave the highest value of ∆K for K = 3 plus a few other suboptimal K values ( Figure S1a) in the analysis of the reduced dataset (one or two representative individuals per population). However, the latter contained clusters with negligible membership ("empty" clusters). Visualization of K = 45 based on the STRUCTURE results (reduced dataset) showed six clusters which are subsets of the clusters in K = 3 ( Figure S1b). STRUCTURE analysis of the whole dataset gave the highest value of ∆K for K = 2 and a suboptimum for K = 5 ( Figure S2). The fastSTRUCTURE results gave a model complexity that maximizes marginal likelihood of 33 (this corresponds to the highest value of ∆K in STRUCTURE). These 33 potential clusters circumscribe mainly the sampling localities/populations with some populations being fused together ( Figure S3). However, both NJ and PCoA analyses based on different distance measures are in correlation with clustering based on STRUCTURE rather than those based on fastSTRUCTURE. The main grouping found in STRUCTURE analyses (K = 2) and PCoA separates the samples into two groups and some admixed individuals ( Figure 2 A deeper look at the clustering patterns in PCoA and STRUCTURE analyses shows that both main groups can be further subdivided. Within the group of the central European highlands and plains (CE), samples from the Havel lowering in Brandenburg form a cluster together with the individuals from Lake Tolk in Schleswig-Holstein, which forms the core CE group. Samples from rivers Aller and Ems in Lower Saxony (CE-AE) appear to be admixed between the core CE group, Danube and Weser. Individuals from rivers Wümme and Eider appear to be related to populations from the Danube region and Weser. The two populations from the river Weser form an individual group with around 1/3 the impact of the British and Rhine populations. Within the second, much bigger group, samples from British rivers form a group as well as the samples from the Baltic and eastern central rivers together with the population from the Theiss river in Hungary (BH). Samples from Danube waters in Austria form a group with more or less impact from the Rhine, Weser and BH. Individuals from waterbodies along the Rhine river are a mixture between the British and the Danube genepools. The same was found for the population from the Botanical Garden of the University in Padua, which should originally be from the Po river. The only population that cannot be assigned to any of the groups is the population from the Danube estuary in Romania because this population shows impacts from Weser, Baltic, CE and Danube genepools.

Principal Coordinate Analyses
Principal coordinate analyses based on the two distance methods (uncorrected p and Dice) gave very similar clustering patterns, with uncorrected p-distances showing a higher total sum of coordinates (Table S1). The first coordinate (uncorrected p: 57%; Dice: 49%) separates the two main groups found in STRUCTURE analyses from each other with the admixed samples in between the two groups. The second coordinate (uncorrected p: 22%; Dice: 25%) separates the two main groups into two subgroups each. (Figure 3). The CE group is separated into the core CE group and the Aller-Ems group (CE-AE). The second group is separated into the BH group and a continuum of samples from British waterbodies, Rhine, Danube, Po and Weser. Among the PCoA based on pairwise F ST distances from hierarchical AMOVAs, those based on groupings according to the STRUCTURE results gave the highest values, and also gave the highest values over all PCoA.

Principal Coordinate Analyses
Principal coordinate analyses based on the two distance methods (uncorrected p and Dice) gave very similar clustering patterns, with uncorrected p-distances showing a higher total sum of coordinates (Table S1). The first coordinate (uncorrected p: 57%; Dice: 49%) separates the two main groups found in STRUCTURE analyses from each other with the admixed samples in between the two groups. The second coordinate (uncorrected p: 22%; Dice: 25%) separates the two main groups into two subgroups each. (Figure 3). The CE group is separated into the core CE group and the Aller-Ems group (CE-AE). The second group is separated into the BH group and a continuum of samples from British waterbodies, Rhine, Danube, Po and Weser. Among the PCoA based on pairwise FST distances from hierarchical AMOVAs, those based on groupings according to the STRUC-TURE results gave the highest values, and also gave the highest values over all PCoA.

AMOVA and Population Statistics
In order to quantify the amount of genetic variation between populations, we have performed AMOVAs. When keeping all sampling sites as separate populations, the analysis showed 97% of the molecular variance occurred between the populations and FST value of 0.97 (Table S1). If populations are assigned according to STRUCTURE results (K = 2), the amount of molecular variance between populations drops to 33% and the FST value to 0.33. To investigate alternative groupings apart from the one based on STRUC-TURE results, we also conducted AMOVAs for groupings based on fastSTRUCTURE results and river systems. Both of these groupings gave higher FST values than the grouping based on STRUCTURE results (Table S1). Average gene diversity over loci in non-hierarchical AMOVA was 0.214; in hierarchical AMOVA based on populations, average gene diversity varied between 0.047 within the commercial samples from Stauden Hameter (population 8) and 0.000 within samples from Potter Heigham (population 21) and Chilley Stream (population 29).

AMOVA and Population Statistics
In order to quantify the amount of genetic variation between populations, we have performed AMOVAs. When keeping all sampling sites as separate populations, the analysis showed 97% of the molecular variance occurred between the populations and F ST value of 0.97 (Table S1). If populations are assigned according to STRUCTURE results (K = 2), the amount of molecular variance between populations drops to 33% and the F ST value to 0.33. To investigate alternative groupings apart from the one based on STRUC-TURE results, we also conducted AMOVAs for groupings based on fastSTRUCTURE results and river systems. Both of these groupings gave higher F ST values than the grouping based on STRUCTURE results (Table S1). Average gene diversity over loci in nonhierarchical AMOVA was 0.214; in hierarchical AMOVA based on populations, average gene diversity varied between 0.047 within the commercial samples from Stauden Hameter (population 8) and 0.000 within samples from Potter Heigham (population 21) and Chilley Stream (population 29).
Nei's H-value (unbiased expected heterozygosity) was estimated with uHe = 0.213 in the overall analysis of all samples together. Analysis of separate populations gave the highest Nei's H-value of H = 0.031 for the commercial samples from Austria (Stauden Hameter population 8) and the lowest value (H = 0.000) for populations from the UK (Potter Heigham, population 21; and Chilley Stream, population 29). Shannon's index was estimated to be I = 0.344 in the overall analysis of all samples together. In the analysis of separate populations, the highest and lowest values were found in the same populations for Nei's H (for details, see Table S2).

Mantel Tests
Mantel tests based on different distance matrices showed between 2.3 and 100% correlation (R 2 ) among the tested pairs of matrices (Table S3). The highest correlations were found between matrices based on uncorrected p-values, Dice distances and binary distances calculated with GenAlEx (r = 0.95-1). The lowest correlations were observed between the matrices based on the genetic data and the matrix containing the geographical data (r = 0.15-0.31), indicating that there is no or only little correlation between the geographic distance and genetic distance of the samples in our dataset. The only pair of datasets where a correlation (r = 0.64) between geographic distance and genetic-based AMOVA distances was observed is the pair of geographic distances and AMOVA distances based on grouping according to river systems.

Discussion
Here, in this study, we examined the genetic diversity of Stratiotes aloides populations from different water systems across Europe. As this species propagates mainly vegetatively, genetic diversity within populations is expected to be low. Due to missing connections (river regulations and lack of flooding) between the water systems, genetic diversity between populations is expected to be high.
Indeed, we did find a high F ST value (0.97) when analyzing the populations separately which shows a high level of genetic differences between the populations and a very low level of genetic differences within the populations, indicating that the populations consist mainly of clones of few genetically different individuals. Considering that the populations propagate vegetatively, and that today, there is no gene flow between the populations via transfer of individuals from one population to another, the relations between the populations could show historical connections between populations. This explains why no, or if only a medium, correlation between geographic distance and genetic distance of the populations is found. The observed correlation between geographic distances and genetic distances based on AMOVA grouped by water bodies has to be viewed with some precaution, as the grouping based on water bodies is, of course, a geography-related grouping and will, therefore, already have a slight bias towards a stronger correlation. Nevertheless, we can still see that there is some correlation between geographic distance and genetic distance when we look at it at the level of waterbodies. The populations from waterbodies of the central European highlands and plains in particular have a genepool which is different from the genepool shared by populations from other regions in Europe. However, it looks like that geneflow between populations has occurred. The fact that water soldier populations from waterbodies along the Rhine seem to be a mixture of genepools from British and Danube genepools might be explained by the historic watercourse of the river Rhine with pervious headwaters of the Danube being directed to the Rhine and a common delta of the Rhine and Thames [23]. A second hypothesis for the connections between populations from British rivers, the Rhine and Italian rivers is long distance dispersal of seeds by migrating water birds (for an example of migration routes of ducks across Europe, see [24]). The connections between the populations from Poland, Baltic countries and the river Theiss in Hungary might also be explained by transfer of plant material by birds [25]. There are not much data available about the dispersal of water soldier fruits by animals, but Efremov et al. [1] and Orsenigo et al. [14] give an overview of the current knowledge of dispersal of Stratiotes and Cook and Urmi-König [2] as well as Forbes [18] mention birds as possible dispersals vectors. Summed up, animals can disperse Stratiotes aloides fruits exo-and endo-zoochorically and if they are migrating over longer distances, seeds and thus genetic information can be transferred between localities. A further point that has to be kept in mind, when investigating relationships among European water soldier populations, is the fact that Stratiotes aloides has a long history as an ornamental plant [1,2]. Unexpected and probably by natural means, difficult to explain relationships between populations could be the result of human-mediated transfer of plant material. As the earliest known fossils of Stratiotes aloides date back around 45 million years [26], the observed groups could be the result of repeated glaciation and deglaciation events in Europe [1]. Summing up, we found the investigated populations of Stratiotes aloides across Europe to form two main groups which can be further subdivided. Roughly, the two groups can be referred to as a central northern Europe-group (CE) and a western-southern-eastern Europe group.
Previous studies of Stratiotes aloides across its distributional range showed a much higher level of genetic diversity within the examined populations [14]. As the sampling regions of the study of Orsenigo et al. [14] are not the same as in our study, the main clustering structures of European populations cannot be fully compared. However, clustering of samples from the Rhine in The Netherlands and the Po in Italy, together with some similarities to populations from the Danube in Bavaria, was observed in both studies. Comparable genetic differences within and between populations of dioecious Hydrocharitaceae were found in Ottelia acuminata where high levels of genetic differences between the investigated populations were found, but little diversity within the populations [27].
All still present-day populations of Stratiotes aloides found in Europe are remnants of much larger and connected populations. For example, in the Danube flood plains around Vienna, Stratiotes aloides was still very common by the mid-19th century, around 100 years later, this species was already mentioned to be rare [28] (and references therein). This example shows that previous large and vital populations became rare and fragmented within the last 150 years.
One possible hypothesis for explaining the differentiation of the samples into two groups could be differences in ploidy level. Orsenigo et al. [14] mention that different ploidy levels (diploid and tetraploid) were observed in Stratiotes aloides. Unfortunately, the material available for our study was not appropriate for chromosome counts and genome size measurements.
Summing up the results and viewing them in light of conservation issues, we can conclude that for conserving the genetic diversity of Stratiotes aloides, it would be preferable to focus on conserving individuals from many different populations all over its distributional range, rather than conserving selected populations with a higher number of individuals per population. For reintroductions, samples from closely located populations, or at least from populations from the same river system, could serve as founder individuals. As sexual reproduction is rare in natural populations, ex situ collections of samples of both sexes might be established to facilitate sexual reproduction and thus, maintain or even slightly increase genetic diversity in Stratiotes aloides. Apart from protecting and conserving Stratiotes aloides as a species, protection of the species as a habitat for fauna species such as the dragonfly Aeshna viridis [29] and the black tern Chlidonias niger [30], which fully or at least mainly depend on Stratiotes [14], is important.

Material
Material was continuously collected between 2012 and 2018 all over Europe wherever populations of Stratiotes were found. Depending on the size of the populations and on the accessibility of the individuals, between 5 and 20 individuals per population were sampled. Wherever possible, individuals from the whole waterbody were collected (e.g., North and South shore, etc.). Short (approx. 5-7 cm) pieces of leaves were collected and immediately dried in silica gel. From several populations, herbarium specimens were collected and deposited in the herbarium of the University of Natural Resources and Life Sciences, Vienna (WHB). Herbarium accession numbers are indicated in the table of accessions (Table 1). In total, we included 345 individuals from 46 populations into the final analysis. As previous studies [14] showed that there is no detectable genetic difference between the two sexes, we did not pay attention to the sex of the collected individuals (for some populations, information about sex is available and can be requested from the authors).

DNA Extraction
DNA was extracted from 20 mg silica gel dried leaf material per individual. The material was ground into a fine powder in 2 mL tubes together with three glass beads in a Tissue-Lyser (Qiagen, Germantown, MD, USA) with 20 s −1 for 5 min. Extraction of DNA was performed via QIAcube (Qiagen) using the DNeasy Plant Mini Kit (Qiagen), mainly according to the manufacturer's protocol. Exceptions were elution of DNA from the columns, which was performed with two steps of 50 µL of elution buffer each. RNA was digested after DNA extraction using 1 µg RNAse A and incubated at 37 • C for 30 min.
Quality control of the DNA extracts was performed photometrically using a NanoDrop 2000 spectrometer. To check RNA digestion, samples were loaded on a 0.8% agarose gel.

AFLP
All DNA extracts that met the quality criteria were adjusted to 100 ng/µL and used for AFLP fingerprinting. Preparation of AFLP samples mainly followed the original protocol [32] with slight modifications.
Restriction of genomic DNA with two restriction enzymes (EcoR I and Mse I) and ligation of double-stranded adaptors to the resulting restricted fragments were performed in one step in a thermal cycler (37 • C for 2 h followed by a hold at 10 • C). Reactions The selective primer pairs (FAM-EcoRI-ACT/MseI-CTA, VIC-EcoRI-ACG/MseI-CTA and NED-EcoRI-ACC/MseI-CTA) were chosen after testing seven different primer combinations in a preliminary test. The selected primer combinations generated clear and not too many bands, thus decreasing the risk of fragments co-migrating by chance, but still with sufficient variability to distinguish the samples.
Reproducibility was checked by repeating ca. 23% of the samples.

Scoring and Phylogenetic Analysis
Sizing and scoring of the data were performed with GeneMarker v2.4.0 (SoftGenetics, State College, PA, USA). After pre-analysis using default settings, sizing profiles of all samples were checked and where necessary, manually corrected. Most of these corrections concerned the 20 bp peak of the size standard. These peaks were often not correctly recognized by the GeneMarker program. High-quality sizing profiles (score > 90) were obtained for all samples. A panel of scorable fragments was established for each primer combination, and fragments between 30 and 600 bp were scored. The relative fluorescent unit (RFU) threshold was set at 40. Automatic scoring was conducted using Local Southern peak call, peak saturation, baseline subtraction, spike removal, pull up correction, and a stutter peak filter of 5% [33]. The results were exported as a presence/absence matrix. The outcome of the automatic scoring was manually checked and corrected for errors. These errors mostly concerned peaks for which shape was atypical. In total, 447 samples corresponding to 345 individuals were scored. From 78 individuals, replicate samples were performed (between two and four replicates per individual). Peak shifts between different analyses dates of the same individuals were used to correct and align all fragment analyses over the whole timespan of the project. These corrections were performed manually and very carefully to avoid artefacts within the dataset. Most of these corrections were small shifts of the majority of peaks by one or two base pairs. For the final analyses, we ended up with 345 individuals, for which high-quality fragment profiles for all three primer combinations could be obtained.
All three primer combinations were combined in a single matrix and analyzed together. Different distance measures were tested for their power to resolve relationships with our dataset. Distance matrixes were calculated in PAUP* v4.0a167 [34] (Nei-Li distance) and SplitsTree v4.15.1 [35] (uncorrected p, Dice, Jaccard and Hamming). Phylogenetic relationships based on previously mentioned distance matrices were reconstructed using SplitsTree to create unrooted NJ dendrograms. To assess the robustness of branches, NJ-bootstrap (NJ-BS) analyses were performed using SplitsTree.
To investigate further significant groupings of the included individuals, we used the programs STRUCTURE v.2.3.4 [39][40][41][42] and fastSTRUCTURE v 1.0 [43]. STRUCTURE was initially run for K = 1-50 with a subset of one or two individuals per population to keep analysis time in a reasonable frame. Based on those results, a second STRUCUTRE analysis with the full dataset (345 individuals) for K = 1-8 was run. We ran STRUCTURE with 10 replicates each and a model based on admixture and independent allelic frequencies, without considering information regarding sampling localities. Each run had 100,000 iterations with 10% additional burn in. The calculation of delta K (∆K) [44] and preparation of the input file for Clumpp were performed with Harvester [45]. Production of a combined file from the ten replicates of the best K was performed using Clumpp v1.1.2 [46] with the Greedy search algorithm. The graphical representation of STRUCTURE results was prepared with Distruct v1.1 [47]. FastSTRUCTURE was ran with the full dataset for K = 1-50. The calculation of ∆K and graphical representation of results were performed with the functions "chooseK.py" and "distruct.py", both implemented in the fastSTRUCTURE package.
Both non-hierarchical and hierarchical analyses of molecular variance (AMOVA) and calculations of population statistics were conducted using Arlequin v3.5.2.2 [38]. The Excel plugin GenAlEx v6.503 [48] was also used for calculating population statistics, AMOVAs, PCoA and Mantel tests. For hierarchical AMOVAs, groups have been defined based on different possible clustering according to populations (sampling locality), river systems (grouped according to Trockner et al. [31]) and STRUCTURE results. Mantel tests [49] were performed based on distance matrices calculated with SplitsTree, pairwise F ST values from AMOVAs, binary distances calculated with GenAlEx and geographic distances (calculated with Geographic Distance Matrix Generator v 1.2.3; [50]). Calculations of Nei's heterozygosity [51], Shannon's information index [52] and percentage of polymorphic fragments was performed with GenAlEx.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/plants10050863/s1, Figure S1a: Delta K values of the reduced dataset (one or two individuals per population), Figure S1b: Visualization of STRUCUTRE results for K = 3 and K = 45 from the reduced dataset, Figure S2: Delta K values of the complete dataset, Figure S3: Visualization of fastSTRUCTURE results for K = 33, Table S1: Results of PCoA and AMOVA, Table S2: Overview of frequency values from population statistics, Table S3: Overview and comparison of Mantel test results. File S1: Data matrix containing AFLP data.

Data Availability Statement:
The data presented in this study are available as supplementary material File S1.