Genetic analysis of two species of Mnomen in the Kalamazoo Watershed reveal panmixia in Z. Aquatica, structure among Z. Palustris, and hybridization in areas of sympatry

Introduction

Wild rice of the genus Zizania, or Mnomen as it is known to indigenous peoples, is both ecologically important to the freshwater ecosystems of Kalamazoo watershed of Michigan and culturally important to native peoples of the Ojibwe Nation. In July of 2010, the Kalamazoo River suffered major damage when a pipeline operated by Enbridge leaked over one million gallons of oil into the surrounding waterways and soil (EPA, 2017). Thirty-five miles of the Kalamazoo River and its tributaries were closed for clean up until June of 2012 (EPA, 2017). Costs to clean up the river and tributaries would eventually exceed $1.21 billion (EPA, 2017). Following clean-up of the oil spill, there was interest in the degree to which the oil spill and subsequent clean-up activities may have impacted Zizania populations along the Kalamazoo River and its tributaries.

Prior to the oil spill and subsequent clean-up, there were three distinct sub-species of Zizania found throughout this system. Zizania palustris is nationally rare, but is commonly found throughout the Great Lakes region of the US. Zizania palustris interior is closely related to Z. palustris, but has recently been recognized as a distinct species. Lastly, Z. aquatica is a riverine or coastal rice species commonly found along large rivers and coastal estuaries of the United States, including a majority of our sample locations in Michigan (Ford-Lloyd, Newbury & Virk, 2001). Wild rice has shown strong evidence for co-adaption in geographically distinct aquascapes that have resulted in co-adaptive gene complexes associated with agronomically important quantitative traits (Jump, Marchnat & Penuelas, 2009).

In 2003, the Anishinaabe people of the Great Lakes region, created an Integrated Natural Resource Management Plan (INRMP) which set goals for the management of cultural and ecologically important resources found on tribal controlled lands (KB-1152-2003). The Wild Rice and Native Plants component of the plan set fourth four key goals relating to native plants, such as Mnomen. These goals centered around the management of the native plants and encouraging increased cultural use. The goals for wild rice include mapping and monitoring populations along natural waterways, as well as monitoring and maintaining adequate adaptive potential (KB-1152-2003). One key knowledge gap noted by the INRMP was the distribution and degree to which all three species co-occurred, or not throughout the Kalamazoo Riverine system.

The goals of this project were to address key knowledge gaps identified by the INRMP. Specifically, we sought to: (1) map the distribution of each of the three species of wild rice along reaches of the Kalamazoo River and tributaries, and specifically focused on areas putatively impacted by the 2010 Enbridge oil spill and subsequent clean-up; (2) assess the existing genetic diversity and structure of wild rice of the Kalamazoo watershed following the oil spill and clean up, and (3) although not specifically expected prior to the start of this study, we also sought to measure rates of hybridization and introgression among extant wild rice along impacted reaches of the Kalamazoo River and its tributaries.

Materials & Methods

Prior to the Enbridge spill in 2010, there were limited to no official records of Mnomen along the main channel of the Kalamazoo River. However, the Anishinaabe regularly harvested Palustris from the Kalamazoo River along the western edge of our study site (Lee Sprig, Anishinaabe Mnomen expert, Pers. comm., 2017–2019 Fig. 1). Following river clean-up, rice beds were also noted east of Albion College—there were anecdotal reports that these beds may have been intentionally transplanted to these locations by some local peoples. Additionally, Aquatica was the dominant and possibly only Mnomen species found along Nottawaseppi Creek prior to the 2010 Enbridge spill. Anishinaabe rice harvesters began to note the occurrence, or an increased occurrence, of both Palustris and Aquatica along Nottawaseppi Creek following clean-up (Lee Sprig, Anishinaabe Mnomen expert, Pers. Comm., 2017–2019). We were directed and accompanied to these various rice beds by local rice harvesters, and tribal members assisted during field surveys.

Prior to the start of field work, 11 locations were identified for sampling based on Traditional Ecological Knowledge of culturally important wild rice beds throughout this system, and access to all sites was made possible by the cooperation of members of the Nottawaseppi Huron Band of the Potawatomi via ceded territories clause of the Keweenaw Bay, Integrated Resource Management Plan of 1936 as amended in 2002–2012. While traveling to the a priori defined sample locations we recorded the location of additional rice beds located, and opportunistically sampled them as well. Using this approach, we sampled 28 beds for a total of 113 rice samples (30 Palustris and 80 Aquatica; Fig. 1).

At each sample location, we collected an ∼2 cm clipping of wild rice leaf tissue and stored it on dry ice until extraction could take place back at the lab. At each sample location we also collected covariate data including coordinates, water temperature, water pH, and putative species identity of rice being collected. In situ species identification was carried out by local indigenous tribal members that are known experts on wild rice identification. Field identification of wild rice was later confirmed using genetic analysis.

Prior to chemical extraction of DNA, all samples were cold shattered by first storing them at −20 °C for 24–72 h, and then mixed with liquid nitrogen and agitated with a Vortex Genie at high speed for 4 h. Samples were then extracted using the EZNA SP Plant DNA kit protocol from Omega BioTek (Omega Bio-tek Inc., Norcross, GA, USA). Zizania samples were screened at 19 polymorphic microsatellite primers that have previously been cross-amplified for Zizania sp. (Richards et al., 2004; Quan et al., 2009; Table 1). Each forward primer was tagged with an M13 universal primer attached to the 5′ end (Schuelke, 2000). We used Kapa 2G Fast HotStart ReadyMix (MilliporeSigma, St. Louis, MO, USA) for PCR assays following manufacturer protocols with the following modifications: each PCR cocktail reaction consisted of the following: 6.25 µl Kapa 2G Fast HotStart ReadyMix, 2.2 µl water, 0.62 µl forward primer, 0.62 µl reverse primer, 0.41 µl M13 primer, 0.7 µl betaine, 1.5 µl DNA (at∼41 ng/µl concentration), for a total volume of 12.3 µl per reaction. Conditions for PCR amplification were as follows: initial denaturation cycle at 95 °C (3 min), then 35 cycles (denaturation, annealing, extension) at 95 °C (15 s)/60 °C (15 s)/72 °C (15 s), followed by one cycle of final extension at 72 °C (3 min), plus an additional 8 cycles of 95 °C (30 s)/53 °C (45 s)/° C 68 °C (45 s), followed by a final extension cycle at 72 °C (3 min) to incorporate the M13 universal primer sequence into the final PCR products. DNA fragments were visualized using the University of Maine DNA Sequencing Facility, Orono ME on an ABI-3730 with Liz-500 size standard.

Multi-locus genotypic data were scored and collated for each individual using GeneMarker 1.97 (SoftGenetics LLC, State College, PA, USA) automated genotypic software (SoftGenetics LLC, State College, PA, USA). The genotypes of each individual were scored manually by three lab technicians and a consensus genotype is used for each individual’s genotype. We initially screened 19 microsatellite markers for use with wild rice in this system. Three markers proved to be monomorphic or did not amplify and were eliminated (Table 1). Examination of allele frequencies for the remaining 16 markers indicated the potential for heteroduplexing in two markers (Dieffenbach & Dveksler, 2003). As a result of heteroduplexing some markers will have higher-than-expected ploidy (Qiu et al., 2001). Conversely, in rare occasions higher than expected ploidy can also result from plant selfing as opposed to biochemistry anachronisms of the PCR reaction. In either case, such markers are generally not suitable for use in standard genotypic analysis of population diversity and structure (Qiu et al., 2001). We also tested for deviations from Hardy-Weinberg Equilibrium, using Program GeneAlEx 6.502 (Peakall & Smouse, 2012), for Linkage Disequilibrium using GenePop on the Web (Rousset, 2008), and for the occurrence of null alleles using Program Microchecker (Oosterhout et al., 2004). These analyses found an additional three markers that were either out of Hardy Weinberg Equilibrium, not at Linkage Disequilibrium or had significant occurrence of null alleles (Table 2). Elimination of these markers left us with 11 microsatellite loci suitable for use in this analysis. After eliminating problematic markers, we also eliminated samples with more than 15% missing data to ensure accurate estimates of gene flow and structure, leaving us with 86 individual samples (29 Z. palustris and 57 Z. aquatica).

We used Program GenAlEx 6.502 to estimate genetic diversity (both expected H_E and observed H_O) within each species and/or subpopulation of rice within each species. We estimated effective population size (N_e) using the Linkage Disequilibrium model (Hill, 1981; Waples, 2006; Waples & Do, 2010), as implemented in NeEstimator V2.1 (Do et al., 2013). We used the method proposed by Jones, Ovenden & Wang (2016) to calculate confidence intervals by jackknifing over individuals rather than loci to control for pseudoreplication associated with the physical linkage among overlapping pairs of loci. We used a Wilcoxon Sign Rank test of allele frequencies to test for a mode shift in allele frequencies commonly associated with a population bottleneck (Soulé, 1987). Additionally, in GenAlEx we calculated allelic richness (AR), and observed and expected heterozygosity H_O and H_E respectively and to test for the Probability of Identity (PI) and the Probability of Identify among siblings (PI_sibs). PI and PI_sibs are generally viewed as measures of the power of the genotypic data set to resolve differences in population structure with a given set of markers.

We used Program Structure 2.2 (Pritchard, Stephens & Donnelly, 2000), to test for hybridization between Z. Palustris and Z. aquatica, and to test for latent population structure within each species (Porras-Hurtado et al., 2013). This analysis works by using Bayesian statistics to group individuals into non-mutually exclusive genetic clusters based on their individual genotypes and the individual allele frequencies within a population (Pritchard, Stephens & Donnelly, 2000). If an individual is assigned equally to two different clusters it would likely be an F1 hybrid of, in this case, Z. palustris and Z. aquatica mating (Hallgren et al., 2003). Thus, structure can identify individuals that are the result of hybridization between two or more species or regional variants. We also mapped the location of putative hybrids on the landscape to gain additional insights into the spatial context or environmental conditions which might promote or inhibit hybridization in this system.

Program Structure was also used to test for population structure within each species. If there is substructure present on the landscape it usually leads to violation of HWE, given that the populations are not panmictic. Geneland (Guillot, Santos & Estoup, 2008), which uses a similar system to Structure but also incorporates spatially explicit data was used to confirm the results of the Structure clustering analysis.

Results

Analysis of genetic diversity with GeneAlEx indicated that both species were highly heterozygous and mostly out of Hardy-Weinberg (Table 2). We also had adequate power to identify population structure and characterize genetic diversity; for Z. aquatica (PISibs <0.0001; H_O = 0.672; H_E = 0.669 AR =13.64) and for Z. palustris (PISibs <0.000081; H_O = 0.636 H_E = 0.697 AR =10.64). The 95% CI for N_e for both rice species encompassed infinity, specifically N_e Z. palustris = 59.5 (95% CI = 22.8-infinity) and N_e Z. aquatica = 663 (95% CI = 114.4-infinity). We similarly detected no significant indication of population genetic bottleneck in Aquatica (Wilcoxon Sign Ran Test P = 0.89697) nor in Palustris (Wilcoxon Sign Ran Test P = 0.68115). Moreover, a plot of the allele frequencies within each species exhibits the predicted L-shaped sigmoidal curve associated with a non-bottlenecked and outbred genetic population.

Analysis with Program Structure to differentiate between Z. aquatica and Z. palustris indicated strongest support for K = 2 genetic clusters corresponding strongly to the field identification of each species. Specifically, 5.8% (five out of 86) of samples were misclassified in the field and identified genetically as different rice species. Of these five misidentified samples, four were identified as Z. aquatica in the field, but genetically identified as Z. palustris, and one was identified as Z. palustris in the field, but genetically identified as Z. aquatica.

When we plotted the cluster assignment of each individual sample using Sigma Plot along with the 95% confidence interval of cluster assignment, we were able to identify five potential hybrid individuals (Fig. 2). Additionally, four out of five (80%) of samples misclassified in the field were located in areas where genetic hybrids were identified or were growing in rice beds with both species present. This suggests that hybridization may have confounded rice expert’s ability to accurately identify some samples.

When we mapped the location of these hybrids, we observe that hybridization tended to occur in areas where both Z. palustris and Z. aquatica are sympatric. Nottawa Creek shows hybrid activity where Z. palustris and Z. aquatica beds are close, as well as sections of the Kalamazoo River upstream from where the oil spill occurred (Fig. 3). In Nottawa Creek, it should be noted that wherever hybrids were detected, both Z. palustris and Z. aquatica were also detected. Downstream of the oil spill in Gull Lake the rice is mostly Z. palustris, with the exception of a couple of small beds of Z. aquatica that occurred in areas where active sediment dredging and new sediment translocation/restoration had occurred following the oil spill (Fig. 4).

Using Program Structure hierarchically, we also tested for substructure within Z. aquatica and Z. palustris. We found equal evidence for K = 1 and K = 3 potential subpopulations for Z. palustris within our study system (Fig. 3), and greatest support for K = 1 of Z. aquatica, indicting a panmictic population. Analysis with Geneland confirmed greatest support for K = 3 genetic subpopulations of Z. palustris (Fig. 4).

The three subpopulations of Z. palustris tend to be genetically diverse and highly heterozygous. Subpopulation three was the most heterozygous and had almost double the genetic diversity of subpopulations one and two. Measures of genetic differentiation as indicated by F_ST values from these three populations show that subpopulation sp1 and sp2 are closely related (F_ST sp1 vs. sp2 = 0.138; sp2 vs. sp3 = 0.141), whereas sp1 vs. sp3 are less genetically linked F_ST = 0.066 (Table 3).

Discussion

Combining Traditional Ecological Knowledge of the location of Mnomen beds along the Kalamazoo River and its tributaries as well as opportunistic encounters, allowed us to sample and map 28 unique rice beds. Prior to this analysis, no official Mnomen beds were recognized along the Kalamazoo River channel (Lee Sprig and Steve Allen, rice experts for the NHBP Tribe, Pers. comm., 2017–2019). Two species were found, Zizania palustris (n = 30) and Zizania aquatica (n = 80), and genotyped them at 11 microsatellite loci. Dredging also likely resulted in the displacement of Z. aquatica to a more southern location than where it was previously suspected of having been present. This indicates that the oil spill may have impacted the presence and distribution of Mnomen in this region.

Of the 86 samples that were genetically analyzed we see a 5.6% misclassification rate, primarily of Z. palustris, and most of these were in areas where hybridization is expected based on the presence of both species. Phenotypic plasticity in both species (Duvall, 1987) may also have led to the likelihood of misclassification in the field. The hybridization of Z. palustris × Z. aquatica leads to a significantly reduced fruit set being produced compared to intraspecific crossing (Duvall, 1987); however, the offspring are fertile. Unnatural sympatry may lead to reduced yields of wild rice, and would be hard to avoid given the predominance of wind dispersal of pollen for these species.

The anthers of Zizania sp. have an independent histology from the flowers and tend to emerge along the lower 1/3 of the 2–3 m long Zizania sp. stem (Zaitchik, LeRoux & Kellogg, 2000). Wind dispersal of pollen is common across all three Mnomen taxa studied (USDA, 2001; MDNR, 2018). However, this region is heavily and densely forested. Consequently, Mnomen tends to have a rather limited dispersal ability (de Wet & Oelke, 1978; Zaitchik, LeRoux & Kellogg, 2000). Therefore, if hybridization among species is to occur, they likely have to be found growing in close proximity to each other. However, hybridization among closely allied lineages is not uncommon among plants when occurring in close proximity to each other (Whitney et al., 2010). Our observations of the occurrence of hybrids only in areas where both species were observed to be growing in adjacent or shared beds further supports the notion that Mnomen has limited dispersal ability. The degree to which hybridization may benefit or negatively impact Mnomen populations of the Kalamazoo warrants further investigation.

Z. aquatica exhibits panmixia and evidence was only found for a single population, whereas for Z. palustris evidence suggests that three populations exist within this stretch of river. Our results suggest that for Z. aquatica gene flow is occurring at a sufficient rate to maintain genetic diversity, though that southern population that was likely displaced by sediment relocation may see differentiation over time if it is not connected to the main populations. The differentiation of Z. palustris into three populations suggest that asexual reproduction may be occurring at a greater rate or there is a barrier to gene flow among these populations, such as geographic isolation which agrees with other studies of this species (Hayes, Stucker & Wandrey, 1989).

While genetic diversity for both species in this study were moderately high it is possible that not enough stands or smaller stands were sampled. Allelic richness was primarily contributed by larger beds for Zizania texana, an endemic and highly geographically restricted wild rice in Texas (Richards et al., 2007). Further sampling may be necessary to determine true populations and genetic diversity as well as to differentiate other potential hybrids. Maintaining or enhancing the presence of larger stands may be an important management goal to sustain local genetic diversity. For example, a study done on Z. latifolia in China suggested that the rivers did not serve as effective corridors for gene flow among populations (Chen et al., 2017). That study showed high inter-population and low intra-population genetic variation was likely the result of inbreeding or clonal breeding as opposed to sexual reproduction necessitated as a result of riverine disturbance resulting in a recent bottleneck (Chen et al., 2017). A similar situation may be evolving in the Kalamazoo River and tis tributaries, such that following the oil spill Mnomen has been displaced to more disjunct and fragmented locations, with limited opportunity for genetic mixing. Genetic monitoring and more high-density surveys are needed. Fortunately, given the cultural importance of the Mnomen to the Anishinaabe People, regular high-resolution surveys of this system occur every year during annual rice harvest, this important reservoir of Traditional Ecological Knowledge about the system is available for resource management agencies to use for future conservation and management planning.

Supplemental Information