Introduction

Hybridization has been associated with increased invasiveness in plants (Ellstrand and Schierenbeck 2000; Schierenbeck and Ellstrand 2009; Hovick and Whitney 2014), which has been attributed to greater phenotypic variation, more extreme trait values, and/or novel phenotypes in hybrids compared to parent taxa (Ellstrand and Schierenbeck 2000; Rieseberg et al. 2007). Overall, naturally occurring invasive hybrids have greater fecundity and are larger than parental taxa (Geiger et al. 2011; Tavalire et al. 2012; Hovick and Whitney 2014). In North American aquatic ecosystems, some of the most aggressive and widespread plant invasions are driven by interspecific or intraspecific hybrids (Galatowitsch et al. 1999; Zedler and Kercher 2004), including hybrid cattail (Typha × glauca; Larkin et al. 2012), invasive cordgrass (Spartina spp.; Ayers et al. 2004; Ainouche et al. 2009), saltcedar (Tamarix spp.; Gaskin and Schaal 2002), and reed canarygrass (Phalaris arundinacea; Lavergne and Molofsky 2007). In these and other cases, one of the parental taxa is itself considered invasive (e.g., Typha angustifolia and Spartina alterniflora), making it important to assess the relative invasiveness of hybrids compared to parental taxa.

Hybrid watermilfoil, Myriophyllum spicatum × Myriophyllum sibiricum (hereafter HWM), in North America is a cross of non-native invasive Eurasian watermilfoil (M. spicatum L., Haloragoraceae; EWM) and native northern watermilfoil (M. sibiricum Kom.; NWM) (Moody and Les 2002, 2007; Zuellig and Thum 2012; Moody et al. 2016). EWM has been present in North America since at least the 1940s and is considered one of the most aggressive—and is among the most highly managed—aquatic invasive species in North America (Aiken et al. 1979; Couch and Nelson 1985; Smith and Barko 1990; Madsen et al. 1991; Bartodziej and Ludlow 1997; Boylen et al. 1999; Kujawa et al. 2017). EWM forms dense canopies at the water’s surface and spreads readily via clonal growth and fragmentation (Aiken et al. 1979; Madsen and Smith 1997)—traits that can cause reductions in native macrophytes (Madsen et al. 1991; Boylen et al. 1999), interfere with aquatic recreation (Aiken et al. 1979; Newroth 1985), and decrease lakeshore property values (Zhang and Boyle 2010). NWM, while superficially similar to EWM, does not exhibit the prolific growth and canopy-formation associated with EWM (Aiken et al. 1979; Aiken and Picard 1980) and is a common, but non-dominant, component of the native flora.

HWM is well documented to pose significant management challenges through the emergence of herbicide tolerance. Specifically, certain strains of HWM are less sensitive than EWM to commonly used herbicides, such as 2,4-d and fluridone (Berger et al. 2012; LaRue et al. 2013b; Netherland and Willey 2017). However, the putative ecological invasiveness of HWM in the field (e.g., the ability to grow fast and large, produce more propagules, and disperse more readily; Van Kleunen et al. 2010; Hovick and Whitney 2014; Moravcová et al. 2015) has largely stemmed from anecdotal observations of high abundance in lakes where the taxa occurs (Moody and Les 2002, 2007). The invasiveness of HWM compared to parental taxa has recently been examined in laboratory experiments, which showed that HWM stems grew significantly faster than parental EWM (LaRue et al. 2013b; Taylor et al. 2017; Thum and McNair 2018). Nonetheless, these relatively short-term, small-scale laboratory studies may have limited transfer to actual lake ecosystems, and we are not aware of any field-based studies that have evaluated the invasiveness of HWM compared to parental taxa. A thorough comparison under natural conditions of key HWM traits compared to those of EWM and NWM is needed to assess the relative invasiveness of this hybrid taxon.

Examination of measures of HWM invasiveness in the field is essential to elucidate ecological impacts and inform management responses (Tavalire et al. 2012). For example, plant size (e.g., growth or biomass)—a key indicator of invasiveness (Hovick and Whitney 2014)—may differ substantially in situ within lakes compared to under laboratory conditions. Surface matting in EWM, a trait that can negatively impact native macrophytes (Madsen et al. 1991; Boylen et al. 1999) and recreational activities (e.g., boating and swimming), is also difficult to examine in the laboratory. Additionally, traits related to reproduction are crucial to examine in the field. Comparing reproductive capacity between hybrid and parental taxa is particularly important for HWM, as hybrid plants are reproductively viable and appear to reproduce sexually in nature (both with other hybrids and through back-crosses with parental taxa) more so than does EWM (Zuellig and Thum 2012; LaRue et al. 2013a; Eltawely et al. 2020; Thum et al. 2020).

Seasonal growth patterns and timing of key life-history events (i.e., phenology) are another important consideration in evaluating the invasiveness of HWM relative to parental taxa. Differences in phenology, such as earlier leaf-out and flowering, can provide invasive species advantages over natives (Godoy et al. 2009; Wolkovich and Cleland 2011; Alexander and Levine 2019). Moreover, the phenology of invasive plants is important for understanding their ecology and management, e.g., when impacts are likely to be most pronounced and when treatments should be applied (Wolkovich and Cleland 2011; Buisson et al. 2017; Wersal and Madsen 2018). The phenology of EWM has been well documented (Patten 1956; Adams and McCracken 1974; Nichols and Shaw 1986; Perkins and Sytsma 1987; Madsen 1993, 1997); however, we know of no studies that have examined the phenology of HWM under field conditions or compared HWM phenology to that of EWM and NWM. Furthermore, while phenological differences have been implicated as drivers of invasiveness in general (Wolkovich and Cleland 2011), this phenomenon has not been well-explored in studies of hybrid invasives relative to parental taxa.

Studies that compare the performance of hybrid taxa to both invasive and native parents are uncommon, but vital for assessing the ecological consequences of hybridization (Hovick and Whitney 2014). Populations of HWM are widespread in lakes throughout the northern United States and Canada (Moody and Les 2007; Borrowman et al. 2014; Grafe et al. 2015; Moody et al. 2016; Eltawely et al. 2020; Thum et al. 2020), making this taxon a pressing concern and priority for management. Improved understanding of the invasive traits and phenology of HWM is needed to determine ecological impacts and guide control efforts. To support this, we tracked populations of HWM, EWM, and NWM year-round in eight Minnesota lakes for 2 years. We measured key indicators of invasiveness and environmental variables that might mediate Myriophyllum performance within each lake. Specifically, we contrasted the relative invasiveness of HWM compared to its parental taxa based on the amount and timing of: (1) flowering, (2) surface cover, and (3) biomass (using stem counts as a proxy).

Methods

Site selection

We used reports from the Minnesota Department of Natural Resources (MNDNR) and previous studies (Moody and Les 2007; Moody et al. 2016) to identify lakes where populations of each Myriophyllum taxa were known to occur within the Minneapolis–St. Paul Metropolitan Area. As the three taxa are visually similar and could be confused in the field, we targeted lakes with only one of the taxa present or dominant (non-target taxa occurred in isolated areas of some study lakes). In spring 2017 we scouted lakes for Myriophyllum beds using an underwater video camera (see Glisson et al. 2020) and marked locations using a handheld GPS unit (Garmin GPSMAP 60Cx; Garmin Ltd.; Olathe, KS, USA). We found sufficiently large populations of Myriophyllum in six lakes in 2017 (Table 1). For each lake, we selected five locations (sites, hereafter) that: (1) were centered within distinct Myriophyllum beds, (2) were distributed throughout the littoral zone, (3) occurred across a range of littoral zone depths, (4) were separated by ≥ 50 m, and (5) were relatively easy to access. In spring of 2018 we repeated this process, incorporating additional locality information provided by Eltawely et al. (2020), and added two additional lakes with five sites each (Table 1). We attempted to identify additional lakes with HWM beds; however, many lakes within the region are treated regularly with herbicides for EWM and HWM control, which limited available study lakes. Indeed, we began sampling one HWM lake, Bald Eagle Lake (Ramsey Co.), in early 2017 but had to remove it from our study following unanticipated herbicide treatment that summer. The remaining lakes included in this study were not subject to whole-lake or large-scale Myriophyllum herbicide treatments during or within 5 years prior to this study (MNDNR, unpublished data). However, when a single HWM plant was identified at a site on Orchard Lake (a NWM lake; Table 1) in 2017, MNDNR performed a spot herbicide treatment at that location. We removed this site from our study and replaced it in 2018 with a new site on Orchard Lake. No other known herbicide treatments were conducted at or within close proximity to sampling sites.

Table 1 Description of study lakes by Myriophyllum taxon, including: taxon (EWM, Eurasian watermilfoil [M. spicatum]; HWM, hybrid watermilfoil [M. spicatum × M. sibiricum]; NWM, northern watermilfoil [M. sibiricum]), lake name, county, geographic coordinates, years sampled, visits, and environmental variables (depth, water temperature [temp.], pH, conductivity [cond.], and dissolved oxygen [DO])

We verified taxon identifications using genetic testing. Myriophyllum samples were collected from all lakes in late August and early September of the year that the lake was first sampled. From each site, we collected multiple apical stems from two plants spaced ≥ 2 m apart. These stems were processed in conjunction with a larger study (Eltawely et al. 2020) and identified to taxon using a genetic assay based on internal transcribed spacer DNA sequences (Thum et al. 2006; Grafe et al. 2015). The genetic assay confirmed the expected morphological identification of each taxon at each site. Further genetic analysis using microsatellite markers (see Eltawely et al. 2020 for details) revealed that genetic diversity within taxa was consistent with previous and concurrent studies (Zuellig and Thum 2012; Eltawely et al. 2020; Thum et al. 2020). That is, EWM had low genetic diversity, with a single genotype represented across study lakes; NWM had high diversity, with eight genotypes across lakes (none shared among lakes); and HWM had intermediate diversity, with three genotypes across lakes (none shared among lakes).

Field sampling

We visited each site every 2–4 weeks during the growing season in 2017 until ice-in (June–November), once in winter 2018 (January or February) and again during the growing season in 2018 until ice-in (May–November). We visited each site on 8–9 occasions in 2017 and 9 occasions in 2018 (Table 1). For each sampling period, we attempted to visit all sites within a five-day window (maximum of eight days between sites for a given sampling period).

Growing-season sampling was performed from a boat. We located sites using a handheld GPS unit and took a new GPS waypoint during each visit to ensure that we were within an acceptable distance of the target location. Upon arrival at each site, we set multiple anchors to maintain our position, centered over the waypoint. We used a 5-m radius from the center of the boat as our sampling unit. Within this sampling unit, we counted all Myriophyllum flower spikes visible above the water’s surface. In 2017, when observing a large number of flower spikes, we used binned ranges (e.g., 75–100) to quantify flower spikes. For analysis, we used the midpoint of the specified range. In 2018, we directly counted all flower spikes except when they exceeded > 500 spikes at a site. In these situations, we counted flower spikes in one quarter of the 5-m-radius circle on opposite corners of the boat, added these two values together, and multiplied this value by two. Sites with > 500 flower spikes had a roughly uniform distribution of spikes throughout the sample space, hence, we considered this subsampling approach adequate for a representative estimate. We estimated Myriophyllum surface cover within the 5-m-radius circle using arsine-square root cover classes (> 0–1, > 1–5, > 5–25, > 25–50, > 50–75, > 75–95, > 95–99, > 99–100; Muir and McCune 1987). All living and apparently rooted Myriophyllum material reaching and/or floating at the surface was included in the measure of surface cover.

As an indicator of biomass, we counted stems at each site by collecting four non-destructive subsamples with a novel underwater video sampling device (briefly described here; see Glisson et al. 2020 for details). We used the sampling device during each visit to a site to record video with an underwater camera (GoPro HERO 5 Black; GoPro Inc.; San Mateo, CA, USA) through the water column at four locations (i.e., subsamples) positioned at fore and aft port and fore and aft starboard, and separated by 2–3 m. In each video, we marked the depth in 0.5-m increments using the GoPro web application on a tablet computer (Apple iPad 2 tablet; Apple Inc.; Cupertino, CA, USA) connected to the camera. We then viewed the videos on a desktop computer and, for each 0.5-m increment, counted the number of stems located between the camera and a 30-cm × 70-cm sampling window positioned 62 cm in front of the camera. We summed the stem counts taken at each 0.5-m interval for each subsample, with the exception of the bottom-most interval (i.e., 0.5 m from the substrate), which could not be collected for all subsamples due to differences in substrate and organic matter on the lake bottom. This resulted in a cumulative stem count metric that, for Myriophyllum in our study lakes, was highly correlated with biomass (Pearson’s correlation coefficient [r] = 0.93; Glisson et al. 2020). To produce a single stem count measurement for each site-visit, we summed the stem counts from the four subsamples and used this value for subsequent statistical analyses.

For winter sampling (January and February 2018), we did not collect flower spike or surface cover data as the lakes were ice-covered. However, we measured stem counts in a similar manner as described above, but from the ice rather than a boat. We located each site and marked locations for two video subsamples located 5-m apart, centered on the site’s GPS waypoint. We then cut a hole through the ice large enough for the sampling device to pass through and recorded and analyzed video samples in the same manner described above (Glisson et al. 2020).

We collected environmental data during each site visit to identify potential differences in growing conditions that could mediate or confound the taxon-specific differences of interest. Specifically, we measured: water depth with our sampling device at two of the four subsample locations around the boat, one on either side of the boat; water temperature (°C), pH, conductivity (µS/cm), and dissolved oxygen (DO; mg/L and percent saturation) at 0.5 m below the water’s surface using a YSI Professional Plus Instrument (YSI Incorporated; Yellow Springs, OH, USA); and Secchi depth using a standard 20-cm black and white Secchi disk.

Statistical analyses

We tested for phenological differences among taxa using counts of flower spikes, proportion of surface cover, and stem counts as response variables. For each response variable, we also plotted the mean value (± 1 SE) from each sampling period for each taxon to visualize general patterns and facilitate comparisons among taxa.

For each response variable, we fit a generalized linear mixed effects model (GLMM) with the glmmTMB package (Brooks et al. 2017) in R statistical software, version 4.0.0 (R Core Team 2020). For surface cover, we converted each cover class to its midpoint, such that each value was a proportion between 0 and 1, and fit a GLMM with binomial errors and logit link. Both flower spikes and stems are count data, so we fit GLMMs for these variables with negative binomial errors and log link. We collected only two stem count samples per site during winter sampling, compared to four samples during the growing season, and we were occasionally unable to measure stems at all four locations around the boat. Hence, we included an offset in the model for stems that indicated the number of samples at each visit (1–4) to account for differences in sampling effort. We were interested in the differences among taxa over time within and between years, so we included the following explanatory variables as fixed effects in models for each response variable: taxon (factor with three levels: EWM, HWM, and NWM), year (factor with two levels: 2017 and 2018), ordinal day of the year (date, hereafter; continuous, starting on January 1: 1–365), a quadratic term for date (date2, continuous), and all two-way interactions among these variables. We included date2 as we expected each response variable to exhibit a curved (i.e., quadratic) relationship over the course of the growing season, with an increase to a peak and subsequent decrease. We scaled date to have a mean of 0 and standard deviation of 1 to improve model convergence and applied this transformation prior to calculating date2. Date2 was highly correlated with water temperature measured at each site visit (r =  −0.90, P < 0.001, for scaled variables), and thus accounts for the relationship between each response variable and water temperatures throughout each year.

To account for repeated measurements from each site, we included a random intercept for site in all models. For each response variable, we also considered a model with random intercepts for lake and for site nested within lake. We compared models with these two random-effects structures using Akaike’s Information Criterion corrected for small sample sizes (AICc) and chose the model with the lower AICc. As a result, we included a random intercept for site in the GLMMs for flower spikes and surface cover, and random intercepts for lake and for site nested within lake in the GLMM for stems. Because each of our response variables included large numbers of zeros, which may indicate zero inflation, we also examined models for each response that accounted for zero-inflation as a function of taxon, year, date, and all additive combinations of these three variables. Once we determined the best random-effects structure, we compared these zero-inflation models with AICc and again selected the model with the lowest AICc. The result was a GLMM for stems that accounted for zero-inflation as a function of taxon and year, and GLMMs for flower heads and surface cover that did not include a zero-inflation term. We then determined the significance of each fixed effect with a Type II Wald chi-square test using the Anova function from the car package (Fox and Weisberg 2019). If taxon was significant in a given model—as a main effect or via an interaction—we calculated the estimated marginal mean value for each taxon (i.e., least-squares mean) and examined overall differences among means using a posteriori comparisons with a Tukey adjustment using the emmeans package (Lenth 2020). We examined model fit and conducted diagnostic tests for models using the DHARMa package (Hartig 2020) and determined that all models were acceptable.

To determine whether there were differences in the environmental variables we measured among the sites selected for each taxon, we fit linear mixed effects (LME) models with each environmental variable as the response and taxon as a predictor using the lme4 package (Bates et al. 2015). We used all environmental data collected during the growing season (May–November). For each environmental variable, we examined a model with a random intercept for site, and a model with random intercepts for lake and for site nested within lake. As above, we compared models with these random-effect structures using AICc, and chose the model with the lower AICc. We then determined whether there were significant differences among taxa in the environmental characteristics of the sites within which they were found. We examined the following environmental variables for each site: water temperature, depth (using mean depth from the two subsamples), Secchi depth (as a proportion of mean depth), pH, conductivity, DO (mg/L), and percent saturation of DO. Some environmental variables were expressed as a proportion or percent (i.e., Secchi depth and percent DO), but their values were not bound between 0 and 1 and their distributions were approximately normal, hence we modeled these variables with an LME model instead of a GLMM. For each environmental variable, we tested for differences in variance among taxa with a Fligner–Killeen test of homogeneity of variances (Conover et al. 1981). If there was significant heterogeneity, we allowed the variance to differ among taxa in the LME model.

Results

Flower spikes

The number of flower spikes significantly differed among taxa and by date2 (Table 2; Fig. 1a). Differences among taxa further differed by year and date2 (significant taxon × year and taxon × date2 interactions, respectively), and the relationship between flower spikes and date2 differed between years (year × date2). Post-hoc analysis indicated that HWM had significantly more flower spikes overall than either EWM (P < 0.001) or NWM (P < 0.001), whereas flower spikes did not differ between EWM and NWM (P = 0.52; Fig. 1a).

Table 2 Model results for fixed-effect explanatory variables from generalized linear mixed effects models (GLMMs) for each response variable (flower spikes, surface cover, and stem counts), including degrees of freedom (df), chi-square values (χ2), and P values
Fig. 1
figure 1

Mean values (± 1 SE) for flower spikes (a), surface cover (b), and stem counts (c) by Myriophyllum taxon. EWM, Eurasian watermilfoil (M. spicatum); HWM, hybrid watermilfoil (M. spicatum × M. sibiricum); NWM, northern watermilfoil (M. sibiricum)

HWM produced more flower spikes overall and maintained consistently more flower spikes throughout the growing season than EWM and NWM (Figs. 1a, 2). Indeed, there were only two sampling periods when a parental taxon had more mean flower spikes than HWM (both in 2017; Fig. 2). The maximum number of flower spikes observed for HWM at a site was 1378 (Otter Lake, August, 2018), compared to 156 for EWM (Cedar Lake, August, 2018) and 109 for NWM (Orchard Lake, August, 2017)—greater than 8- and 12-fold differences, respectively (Fig. 2). Flower spikes of HWM were observed earlier than those of EWM and NWM in 2017 and earlier than NWM in 2018 (Fig. 2). HWM also appeared to produce flower spikes faster, having greater numbers than either parental taxon on the first visit when we observed flower spikes each year. For all taxa, flower spikes appeared by early to late June, but peak flowering times differed—EWM peaked in August through mid-September while NWM peaked in July through the end of August (Fig. 2). The peak in flower spikes for HWM was higher than either parental taxon in both years. In 2017, two peaks were observed for HWM (Fig. 2)—one in late June and one in mid-September; this was primarily driven by sites from Otter Lake, which had both early- and late-season flowering peaks. In 2018, HWM flowering peaked mid-August through mid-September. Lastly, we observed flower spikes of HWM and EWM later in the growing season than those of NWM; this was true in both years but more pronounced in 2018 (Fig. 2).

Fig. 2
figure 2

Flower spike phenology by Myriophyllum taxon in 2017 and 2018. Points and error bars are means ± 1 SE for each sampling period. Note log10-transformed y-axis; all points at the bottom of the plot have a value of 0 (undefined on the log scale). EWM, Eurasian watermilfoil (M. spicatum); HWM, hybrid watermilfoil (M. spicatum × M. sibiricum); NWM, northern watermilfoil (M. sibiricum)

Surface cover

Surface cover significantly differed among taxa and by date2 (Table 2; Fig. 1b). The relationship between surface cover and date significantly differed among taxa (significant taxon × date interaction), and there was modest evidence that the relationship between surface cover and date2 differed between years (year × date2; Table 2). Post-hoc analysis among taxa indicated that HWM had significantly greater surface cover overall than NWM (P = 0.003), but not EWM (P = 0.16); surface cover for EWM and NWM marginally differed (P = 0.063; Fig. 1b).

Surface cover reached greater annual peaks and was greater overall for HWM than for both parental taxa (Figs. 1b, 3). HWM was recorded at the highest observed surface cover (0.85) on 10 occasions, whereas EWM and NWM never reached this level at any site over the 2 years of observations; EWM and NWM reached maximum surface covers of 0.65 and 0.375, respectively. Both HWM and EWM reached the surface earlier than NWM in both years. We observed differences in the timing of peak surface cover among taxa and between years. Surface cover of HWM peaked later in the season than EWM, with peaks in early September in 2017 and mid-July through mid-August in 2018 (Fig. 3). In 2017, HWM appeared to also have an early-season peak in early July; the drop in mid-July was driven by low surface cover at all sites on Thomas Lake, which was likely correlated with flowering phenology, as flower spikes drop below the surface once pollinated and slough from the plant when fruits mature (Patten 1956; Aiken et al. 1979). Surface cover at all sites on Thomas Lake rebounded after this sampling period. Surface cover for EWM peaked mid-June through mid-July in both years, followed by a sharp decline (Fig. 3). In both years, HWM maintained greater surface cover much longer into the fall, with substantially more surface cover than both parental taxa from mid-August through the end of the growing season. NWM had peaks in August in both years and had substantially less surface cover overall than the other two taxa (Figs. 1b, 3).

Fig. 3
figure 3

Surface cover phenology by Myriophyllum taxon in 2017 and 2018. Points and error bars are means ± 1 SE for each sampling period. EWM, Eurasian watermilfoil (M. spicatum); HWM, hybrid watermilfoil (M. spicatum × M. sibiricum), and NWM, northern watermilfoil (M. sibiricum)

Stems

Stem counts significantly differed between years, and differences among taxa differed between years (significant taxon × year interaction; Table 2). Stem counts significantly differed by date2 and this relationship also differed between years (year × date2; Table 2)—the quadratic relationship was concave upwards in 2017 and nearly flat in 2018 (Fig. 4). Post-hoc analysis indicated that stem counts overall did not significantly differ among taxa (EWM–HWM: P = 0.55; EWM–NWM: P = 0.66; HWM–NWM: P = 0.96; Fig. 1c).

Fig. 4
figure 4

Stem count phenology by Myriophyllum taxon in 2017 and 2018. Points and error bars are means ± 1 SE for each sampling period. EWM, Eurasian watermilfoil (M. spicatum); HWM, hybrid watermilfoil (M. spicatum × M. sibiricum); NWM, northern watermilfoil (M. sibiricum)

While there were not statistically significant differences in stem counts among taxa, HWM did produce more stems than either parental taxon (Fig. 1c) and the maximum stem count for HWM was 379, well exceeding that of NWM (255) or EWM (248). There were more stems overall in 2017 than 2018 and the pattern of stems over time differed between years (Fig. 4). In 2017, there was a relatively large number of stems early in the season for all three taxa, particularly for HWM and EWM. For both EWM and HWM, stem counts decreased sharply, leveled off, and then increased through the end of the growing season in 2017. Stem counts for NWM remained steady and then increased from mid-July through the end of the 2017 growing season. In general, stem counts for all three taxa peaked near or at the end of the 2017 growing season in October and November. The patterns in 2018 were less clear. We observed healthy stems for all three taxa under the ice in the winter of 2018. For HWM, there was a general increase in stems throughout the growing season in 2018, with a peak at the end of the season in November, similar to the late-season peak in 2017 (Fig. 4). In 2018, stems for EWM peaked in late June and those for NWM peaked in early September.

Environmental variables

We did not find any significant differences in the environmental variables we examined among the sites occupied by the three taxa, including water temperature, depth, Secchi depth, pH, conductivity, and DO (mg/L and percent saturation; Table 1). However, pH, conductivity, and DO were generally lower in HWM lakes, Thomas Lake in particular (Table 1). Site depth was also shallower on average for HWM lakes, but again, this was largely driven by Thomas Lake. Sites with NWM had the greatest mean depth, with two of the three NWM lakes having the greatest mean site depths among all study lakes.

Discussion

We found field-based evidence that hybridization between invasive EWM and native NWM increases invasiveness with respect to ecologically relevant traits and phenology. Namely, compared to both parental taxa, HWM showed substantially greater reproductive potential and surface cover throughout the season. This reflects differences in plant architecture among the three taxa; HWM did not consistently produce more biomass than its parental taxa, but produced more lateral stems and flower spikes on the water’s surface. To our knowledge, this is the first study to demonstrate increased invasiveness of HWM under natural conditions, and the first to systematically document phenology in HWM. Our findings show key phenological differences among taxa; specifically, HWM produced more flower spikes earlier and maintained consistently more flower spikes throughout the growing season, and HWM had higher peaks in surface cover and maintained higher surface cover later in the growing season. Importantly, the differences we observed among taxa are unlikely to have been an artifact of environmental differences among study lakes, as site depth and water conditions were similar. Other environmental factors not assessed in this study, such as sediment texture and nutrient availability, affect growth and biomass and should be explored in future studies. Nonetheless, these taxa appear to occupy generally similar lakes across Minnesota (Eltawely et al. 2020).

Our study reinforces evidence for increased invasiveness in HWM derived from laboratory experiments (LaRue et al. 2013a; Taylor et al. 2017; Thum and McNair 2018) and adds to a larger and growing body of evidence showing greater invasiveness in hybrid plants (Schierenbeck and Ellstrand 2009; Hovick and Whitney 2014) and substantial phenological differences between invasive and native taxa (Wolkovich and Cleland 2011). Our findings particularly provide support for one component of the hybridization-invasion hypothesis examined by Hovick and Whitney (2014): that hybridization is associated with increased fecundity in invasive taxa.

Our results demonstrate substantial differences among taxa in traits associated with invasiveness; however, traits of hybrid taxa can vary between generations and lineages (Whitney et al. 2006; Hovick and Whitney 2014). Phenotypic variation among HWM lineages has been documented (Berger et al. 2012; Netherland and Willey 2017; Taylor et al. 2017) and there was underlying genotypic variation in HWM and NWM—both within and among study lakes. While our study did not account for intraspecific variability, laboratory studies that have examined multiple HWM genotypes suggest that increases in HWM growth rate and biomass are general phenomena (LaRue et al. 2013b; Taylor et al. 2017). Hence, we believe our findings are likely general to HWM and not specific to the genotypes in our study. Nonetheless, it would be beneficial in future field studies to investigate intraspecific variation in phenology and invasiveness in these taxa, which would be facilitated by recent information on the diversity and biogeography of Myriophyllum genotypes (Eltawely et al. 2020; Thum et al. 2020).

The ability of HWM to produce substantially more flower spikes than EWM and NWM, produce more flower spikes earlier, and maintain consistently more flower spikes throughout the growing season has implications for further spread of HWM. While we did not evaluate reproductive viability, HWM is known to produce viable pollen and seeds, and its high genetic diversity across the landscape is consistent with sexual reproduction being a major component of its expansion (Zuellig and Thum 2012; LaRue et al. 2013a). HWM seeds have been found to have relatively high germination rates (> 50%; LaRue et al. 2013a; Thum and McNair 2018), exceeding those of NWM and being similar to germination rates of EWM (Thum and McNair 2018). While there are relatively few studies on the reproductive biology of HWM, studies of EWM offer additional insights. Pollination rates for EWM (measured by seed set) are high (> 70%; Wani and Arshid 2013) and seed viability for EWM (germination rate) is relatively high under favorable conditions (Patten 1955; Madsen and Boylen 1989; Hartleb et al. 1993; Standifer and Madsen 1997; Wani and Arshid 2013). Total seed set for EWM may be as high as 4 million seeds per hectare in dense populations (Aiken et al. 1979). EWM seeds can retain their viability for several years (Guppy 1897; Patten 1955) and are highly resistant to desiccation (Patten 1955; Standifer and Madsen 1997), providing opportunities for greater spread. These reproductive traits require further investigation in HWM. However, given that each EWM flower spike contains ~ 30 pistillate flowers and given HWM’s and EWM’s high pollination and germination rates (Madsen and Boylen 1989; LaRue et al. 2013a; Wani and Arshid 2013), the quantity of HWM flower spikes we observed suggests a high capacity for HWM to reproduce sexually and spread via seeds.

The ecological consequences of sexual reproduction in HWM are important to consider for both the spread of HWM and impacts on parental taxa. HWM may cross within and among HWM populations or with either parental taxon found in the same waterbody, as demonstrated by the frequency and extent of later-generation hybrids (e.g., F2 and backcrosses; LaRue et al. 2013a). Introgression with both NWM and EWM has been reported for HWM (Zuellig and Thum 2012; LaRue et al. 2013a). The increased ability of HWM to sexually reproduce increases the probability of introgression with native NWM, which could result in displacement of NWM in lakes where both taxa occur. The fact that many Minnesota lakes appear to now contain only HWM (Eltawely et al. 2020) suggests that such displacement via introgression may have already occurred in some locations. Moreover, continued HWM reproduction could produce later-generation hybrids with greater invasiveness than F1 individuals (Whitney et al. 2006; Hovick and Whitney 2014). The extent of crosses within and among taxa, and the genetic makeup of Myriophyllum populations, is an area of ongoing research. Continued genetic monitoring over time should help elucidate the long-term effects of increased reproductive capacity on HWM and parental taxa.

The greater surface cover in HWM throughout the year likely facilitates greater vegetative spread compared to EWM and NWM and may increase effects on native macrophytes. Substantial surface cover increases the likelihood of Myriophyllum fragmentation by watercraft, wave action, or other disturbances, with resulting fragments able to produce new plants (Kimbel 1982; Evans et al. 2011; Thum et al. 2017; Heidbüchel and Hussner 2019). Indeed, fragmentation is considered a primary mode of spread for EWM (Smith and Barko 1990; Madsen and Smith 1997). The relative production of HWM fragments compared to EWM and NWM, as well as the relative capacity of HWM fragments to remain viable and propagate in the field, requires further investigation. We did not evaluate the capacity of HWM to spread via fragmentation, but given its ability to propagate from fragments following herbicide application (Thum et al. 2017), and the substantial quantity of HWM plants we observed at the surface throughout the growing season, HWM appears to have high potential to spread in this manner.

Increased surface cover throughout the growing season could also result in greater negative impacts of HWM on native macrophytes through shading, as has been observed in EWM (Madsen et al. 1991). While HWM did not have consistently greater surface cover early in the season compared to EWM, it did achieve greater surface cover and retain it longer through the year than either EWM or NWM. This greater phenological niche breadth than its parental taxa could allow HWM to access and appropriate resources for longer throughout the year (i.e., the niche breadth hypothesis; Wolkovich and Cleland 2011). EWM has substantial phenological overlap with native macrophytes, increasing potential for competitive displacement of native species by EWM (Verhoeven et al. 2020). Hence, the even broader phenological niche for HWM could exacerbate such effects. In particular, greater phenological niche breadth and shading ability, combined with introgression through increased reproductive ability, may lead to declines in NWM in lakes where both taxa occur.

We did not find strong evidence that biomass, as measured by stem counts, differed among taxa. Interestingly, NWM produced consistently greater biomass than either of the invasive taxa in 2017 and greater biomass than EWM across both years. These findings may seem to contrast with experimental studies demonstrating greater HWM growth compared to EWM and NWM (LaRue et al. 2013b; Taylor et al. 2017; Thum and McNair 2018); however, plant growth in a controlled setting, particularly on single unrooted stem fragments, is likely to differ from growth under natural conditions. Moreover, biomass gained in time-limited experiments with single stems may be poor proxies for biomass of whole plants in situ within and across growing seasons. For example, the biomass associated with stem counts in this study reached 204 g for a single subsample (2–4 subsamples were collected during each site visit), with a mean of 31 g (Glisson et al. 2020), whereas the greatest mean biomass among Myriophyllum taxa documented by Thum and McNair (2018) was < 3 g. Other field studies have demonstrated similarly high levels of Myriophyllum biomass (Adams and McCracken 1974; Madsen 1997). Furthermore, mesocosm experiments have shown that biomass is generally similar among the three Myriophyllum taxa. Specifically, HWM from different populations produced similar or lower biomass than EWM (Berger et al. 2012) and EWM produced similar biomass to NWM (Valley and Newman 1998).

Instead of producing more biomass as a whole, HWM (and EWM) may simply allocate greater biomass to lateral branching at the water’s surface (rather than to main stems in the water column) than NWM—a difference in plant architecture rather than net biomass. Indeed, this variation in architecture is regarded as a key difference between EWM and NWM (Aiken et al. 1979; Aiken 1981). Such structural modification has also been demonstrated in an aggressive strain of the invasive aquatic plant, Hydrilla verticillata, which may also have resulted from hybridization (Benoit et al. 2019); monoecious H. verticillata exhibits greater lateral shoot and tuber formation, despite having similar biomass, than the less invasive dioecious strain (Steward and Van 1987; Van 1989; True-Meadows et al. 2016). While we did not quantify lateral branching, HWM has been shown to branch more than NWM and EWM (Thum and McNair 2018) and we observed the highest surface cover in HWM. Investment in greater light-capturing ability via increased plant height—or in the case of HWM, greater surface cover—is a key trait associated with invasiveness (Pyšek and Richardson 2007; Moravcová et al. 2015) and greater competitive ability in general (Gaudet and Keddy 1988). Increased lateral branching and associated surface cover also facilitates the production of more flower spikes, as these structures occur exclusively on branches at the water’s surface, further increasing the spread potential of HWM.

While prior mesocosm studies and consideration of plant architecture support our findings regarding stem count differences among taxa, our measurements may nonetheless have been affected by our sampling strategy. Myriophyllum beds can appear continuous at the water’s surface, but individual plants are often ≥ 0.5 m apart at the substrate. We did not target individual plants with our sampling device and often measured zero stems for a given sample (hence the zero-inflated data distribution) despite the proximity of multi-stemmed plants. More targeted approaches, with individual plants sampled or randomly selected plants harvested on each visit, may better ascertain biomass differences among taxa in situ.

EWM is considered highly invasive and is a focus of management throughout the U.S. Our findings suggest that HWM may be more invasive in some respects. Furthermore, management of HWM is made more challenging with the hybridization-driven emergence of herbicide tolerance in this taxon (Berger et al. 2012; LaRue et al. 2013b; Parks et al. 2016; Netherland and Willey 2017; Nault et al. 2018). Our description of HWM phenology could help refine management tactics. For example, mechanical harvesting could be used early in the growing season to reduce HWM flower production, or be implemented continuously throughout infested areas to keep plants from matting on the surface (Crowell et al. 1994). However, mechanical harvesting can spread fragments and propagules within a lake (Hussner et al. 2017), an unintended consequence that needs to be accounted for. Herbicide treatments later in the season—when native plants have already senesced but HWM remains active—may also be an option in some situations. Lastly, as certain HWM lineages are less sensitive to herbicides (Berger et al. 2012; Parks et al. 2016; Netherland and Willey 2017; Nault et al. 2018), knowledge of not only taxon (HWM versus EWM) but also HWM genotype could contribute to improved management.

Management of HWM is driven by multiple considerations, including the taxon’s interference with recreation and other human uses. However, management of HWM to reduce recreational and ecological impacts should be considered in light of these impacts and potential non-target impacts of treatment. Recent reports indicate that treatment of both EWM and HWM can result in declines in native macrophytes without desired reductions in target species (Nault et al. 2014, 2018; Kujawa et al. 2017; Mikulyuk et al. 2020). Long-term treatment of invasive Myriophyllum populations may also select for more herbicide-resistant strains of HWM (LaRue et al. 2013b; Parks et al. 2016). The impacts of HWM invasion could be exacerbated if herbicide-resistant strains of HWM also share the invasive traits of HWM observed in this study, potentially resulting in parallel increases in spread and impacts to native macrophytes. While the ecological impacts and effects of management on parental EWM are well documented, these remain critical areas for further HWM research to understand the consequences of this taxon’s increased invasiveness.