Population structure of three invasive congeneric teasel (Dipsacus) species

Abstract

Three species of the Old World genus Dipsacus L. are considered invasive in the Americas, yet they may differ in how they spread and reproduce and in their genetic diversity. Differences in invasion method may suggest that different management techniques are needed for each species. We performed genetic analyses on 572 plants in 69 populations from the United States, Argentina, and Eurasia with the goals of analyzing taxonomy, diversity, mode of reproduction, population structure, and founder effect of each of these species’ invasions, as well as looking for evidence of recent or ongoing hybridization. We found Indian teasel [Dipsacus sativus (L.) Honck.] to be lowest in diversity and possibly reliant on self-pollination more than the other species, Fuller’s teasel (Dipsacus fullonum L.) and cutleaf teasel (Dipsacus laciniatus L.). We found no evidence of hybridization within the invasions and no support for D. sativus as a subspecies of D. fullonum. The closest genetic matches of D. fullonum from the United States to the native range were with Hungary and Spain, while the closest match for D. fullonum between Argentina and the native range was with Spain. Dipsacus laciniatus from the United States most closely matched with samples from Russia. Population structure information regarding these three weedy Dipsacus species can help us understand their invasive processes as well as give insight into their management and the development of a biological control program.

Management Implications

Teasels (Dipsacus species) can form large monocultures, cause loss of riparian area integrity, and occupy habitats important to sensitive or threatened plant species. There are three nonnative teasels in the United States. Different weed species, even within the same genus, can invade differently and may require different control methods. To better understand each of the teasel invasions, we used genetic analysis and found that Dipsacus sativus (Indian teasel) primarily relies on self-pollination while Dipsacus fullonum (Fuller’s teasel) and Dipsacus laciniatus (cutleaf teasel) primarily outcross. We found no evidence of hybridization between species, although this has been suggested from morphological analyses. We also found the closest genetic matches between invasions and the native range, which informs searches for biological control agents.

Introduction

Invasive plant species, even congenerics, can vary in how they spread and persist, and thus may require different management strategies (Mortensen et al. 2000). Differences in reproductive mode, plasticity, phenology, trophic interactions, and abiotic and biotic resistance and tolerance partially drive invasiveness (e.g., Gerlach and Rice 2003; Hao et al., 2017). Understanding these traits can inform effective control methods for existing and new populations (Byers et al. 2002). Additionally, any interspecific hybridization, especially if novel or between native and nonnative congeners, may create individuals that invade differently from parental species (e.g., Grosholz 2010; Larkin et al. 2012; Mayonde et al. 2016) and may also require control and management methods that differ from those used on the parental species (Gross and Rieseberg 2005; Moody et al. 2008; Williams et al. 2014).

There are multiple teasel species in the genus Dipsacus (Caprifoliaceae family; formerly in the family Dipsacaceae) listed as invasive in North America; Fuller’s teasel (Dipsacus fullonum L.), cutleaf teasel (Dipsacus laciniatus L.), and Indian teasel [Dipsacus sativus (L.) Honck.], and it should not be assumed that they all invade in the same manner. As a group they are widespread across the United States, only absent from Alaska, Hawai’i, North Dakota, Louisiana, and the extreme southeast (South Carolina, Georgia, and Florida), but are less common in the Great Plains and desert regions; in Canada, they are present mostly in the southeastern and southwestern provinces (iNaturalist n.d.; USDA-NRCS 2023). The species have different invasive ranges in North America (Figure 1): Dipsacus fullonum is the most widespread, while D. laciniatus occurs mostly in the eastern half of North America, and D. sativus is most numerous in California and the northeastern United States. No Dipsacus species are native to North America, and D. fullonum and D. laciniatus are listed as invasive by 16 states (Rector et al. 2006), where they outcompete many native species (Werner 1975). They are cited as having negative ecological effects, such as development of large monocultures (Weber 2003), loss of riparian area integrity (Ringold et al. 2008), and occupation of habitats important to sensitive or threatened plant species (Snyder and Kaufman 2004), and are listed as invasive in four U.S. national parks (USDI-NPS 2003). Teasel establishment and spread are common on disturbed sites but may also occur in established vegetation (Solecki 1993) and natural areas (Hilty 2009). Another taxon, Dipsacus sylvestris Hudson, is considered a synonym of D. fullonum (Ferguson and Brizicky 1965). Dipsacus fullonum is also invasive in the Pampean region of Argentina (López-Lanús 2016), where it is considered an alternative host for sunflower chronic mottle virus (Giolitti et al. 2009). The center of origin of the invasive teasels appears to be southern Europe (Verlaque 1985), although most are also found in temperate Asia and northern Africa (Weber 2003).

Figure 1. Distribution of invasive Dipsacus species in North America (USDA-NRCS 2023).

Teasel is used in bird seed mixes (Topham 1968) and in flower arrangements for cemeteries (Bentivegna 2006; Bentivegna and Smeda 2011a, 2011b), with both activities likely being sources of teasel invasion. Dispersal along roadways, waterways, and urban expansion is also important to its spread (Skultety and Matthews 2017; Werner 1975). Dipsacus sativus has historically been selected for receptacle bracts that are stiff and recurved to effectively raise the nap on cloth and wool. It was used as such since Roman times and was a popular crop in England in the 14th century (Topham 1968) until more recently, when cultivation moved to France, Spain, and Italy. Dipsacus sativus may have been introduced to North America as early as the 1700s (Donaldson and Rafferty 2002), with reports of cultivation in New York (1840) and Oregon (1907), USA (Dallimore 1912); and it was still under cultivation in California in the mid-20th century (Rector et al., 2006; Stoner 1951). Dipsacus fullonum and D. laciniatus do not have receptacles suitable for raising nap on cloth, but D. fullonum is commonly named Fuller’s teasel (a fuller is a person who works with cloth). This confusion of common names is likely due to D. sativus once being listed as a subspecies of D. fullonum and sharing the common name of Fuller’s teasel. The teasel species that are not optimal for textile processing may have been introduced accidentally with D. sativus (NISC 2023).

Dipsacus fullonum and D. laciniatus are for the most part outcrossing and protandrous, are not known to propagate new ramets from vegetative material, but can self-pollinate at low rates (Bentivegna and Smeda 2011b; Gucker 2009; Verlaque 1985; Werner 1975). They are considered biennials but may stay as rosettes for more than 1 yr and are thus at times considered monocarpic perennials (Gross 1984). There are reports of hybrids between the three invasive species, but plants having intermediate morphological characteristics are found only rarely, and no hybrids have been officially named (Werner 1975). All three species have a diploid chromosome number of 2n = 18 (Temsch and Greilhuber 2010). Control of teasels is currently limited to mowing, herbicide applications, and revegetation (Bentivegna and Smeda 2012; Daddario et al. 2021; Dudley et al. 2009); an investigation into biological control (Rector et al. 2006) was initiated but is currently not progressing.

Our goals are to use molecular markers to investigate the diversity, population structure, and founder effect of each of these species’ invasions, to determine dominant mode of reproduction, and to look for evidence of recent or ongoing hybridization. We also investigate the taxonomic hypothesis that D. sativus is a subspecies or variety of D. fullonum and compare invasive and native genotypes of the three taxa to elucidate invasive species origins.

Materials and Methods

We collected young, disease-free leaves from 572 plants in 69 populations from the United States (n = 298), Argentina (n = 54), and Eurasia (n = 220, primarily Europe) (Figure 2; Table 1; Supplementary Data File, Population data tab) with a range of 7 to 10 (mean of 8.2) plants per population. Some additional collections were just one plant per location, and these were not included in any population-level analyses. We haphazardly sampled plants at least 5 m apart in each population and stored leaves in silica desiccant at ambient temperature. When collecting, we identified plants to species using these key features (Illinois Wildflowers 2023; Jepson Flora Project 2023):

1. Pinnatifid leaves ………… D. laciniatus

1. Entire or toothed leaves … 2

2. Erect or upcurved involucre bracts; receptacle bracts ± flexible, ending in straight spine …………… D. fullonum

2. Spreading or reflexed involucre bracts; receptacle bracts are very stiff, ending in recurved spine … D. sativus

Figure 2. Plant collection locations for Dipsacus fullonum (blue), Dipsacus sativus (green), and Dipsacus laciniatus (red) from (A) the United States, (B) Argentina, and (C) Eurasia. Population labels are noted next to symbols. The one red/blue symbol in A represents a population that contained both D. fullonum and D. sativus.

Table 1. Dipsacus collection and genetic analysis information. ^a

^a G = number of unique amplified fragment length polymorphism (AFLP) genotypes; G/N = number of unique genotypes divided by number of plants sampled; PLP = percentage of loci that are polymorphic at the >5% level; Φ_PT = percentage of molecular variance among populations. Lowercase letters indicate significant differences among populations of the three species in the United States (first letter), and within D. fullonum and D. laciniatus (lowercase letters in parentheses) between United States (invaded range) and Argentina (invaded range) or Eurasia (native range).

We extracted genomic DNA from approximately 20 mg of leaf material using a modified CTAB method (Hillis et al. 1996). The amplified fragment length polymorphism (AFLP) method followed Vos et al. (1995) with modifications as in Gaskin and Kazmer (2009). All 15 selective primer combinations of MseI + CAA, CAC, CAT, CTA, or CTC and EcoRI + AAG, ACC, or ACT were prescreened for PCR product quality and number of variable loci using eight samples, and the two most polymorphic primer pairs were chosen (viz., MseI + CAC/EcoRI + ACT and MseI + CAT/EcoRI + ACT). We omitted AFLP data from any plants that did not produce a typical electropherogram pattern (i.e., noise >20 relative fluorescence units [rfu] or failure to produce peaks). We made final allele calls for loci manually with ABI GeneMapper (ThermoFisher Scientific, Waltham, MA, USA) at >50 rfu; bin width of 1 bp.

We performed DNA sequencing of the nuclear ribosomal internal transcribed spacer (ITS) region as in Gaskin et al. (2020) for 11 plants (D. fullonum, n = 6; D. laciniatus, n = 2; D. sativus, n = 3) using the forward and reverse primers ITS 1 (5′-TCCGTAGGTGAACCTGCGG-3′) and ITS 4 (5′-TCCTCCGCTTATTGATATGC-3′) from White et al. (1990). We aligned sequences in MEGA X (Kumar et al. 2018), and a haplotype network was constructed manually. DNA sequences are listed in the Supplementary Data File, Sequence data tab.

We calculated Dice pairwise similarities to assess genetic similarity within populations of each Dipsacus species in both the native and introduced range. Genetic similarity (Dice: 2a/(2a + b + c), where a is the number of bands present in both samples and b and c are the number of bands present in only one or the other sample, respectively) between genotypes was calculated using the DIS/SIMILARITY module of NTSYS-pc v. 2.1 software (Rohlf 1992). To estimate AFLP PCR error rate, we performed repeats of 48 plants (8.4% of the total 572 plants) starting with CTAB-extracted material, scored them blindly, and calculated the number and percentage of mismatches between the original and repeat AFLP data sets. We counted the number of genotypes (G) in a population manually in a spreadsheet of Dice similarity values. Under the assumption that an increase in identical genotypes in a population indicates less outcrossing and more self-pollination, we compared the mean proportion of unique genotypes detected in populations of each of the three species and in the invaded versus the native range. Data were analyzed in R v. 4.3.1 (R Core Team 2023). We used binomial generalized linear models compared with type II ANOVAs (function Anova in the car package) followed by post hoc Tukey tests (function emmeans in the emmeans package) to assess mean differences in G/N (number of unique genotypes out of number of plants sampled) among populations. We examined differences among U.S. populations for which we had at least seven samples of the three species, and between native and invaded regions for each of the two species for which both native and invaded range data were available. We calculated proportion of loci that are polymorphic (PLP) at the ≥5% level manually in a spreadsheet.

To visualize clustering of AFLP genotypes, we performed principal coordinates analysis (PCoA) using Dice values and the DCENTER and EIGEN modules of NTSYS-pc for all three species combined and for each species separately. To determine the number of genetic clusters (K) represented in the genotypes, we performed population clustering and assignment tests using the software STRUCTURE v. 2.3.3 (Falush et al. 2003, 2007; Pritchard et al. 2000). Binary AFLP data were diploidized (i.e., no peak at a locus was scored as 0/0; peak at a locus was scored as 1/unknown, because AFLPs are dominant data and thus ambiguous if presence = 1/1 or 1/0 when coding for codominant data input; see Falush et al. 2007), no population or geographic location information was included, admixture was assumed as possible, allelic frequencies were considered to be independent, and a 50,000-run burn-in (α stabilized at approximately 1,000 runs) and 100,000-run length were used. We tested for number of genetic clusters (K = 1 to 10) with 10 repetitions for each value of K. Selection of K from these output data was done with the criterion ΔK suggested by Evanno et al. (2005), and results were visualized in the software STRUCTURE HARVESTER web v. 0.6.92 (Earl and vonHoldt 2012).

To analyze population structuring we performed distance-based analysis of molecular variance and resulting genetic differentiation (Φ_PT) on the binary AFLP data, using the GenAlEx add-in for Excel (Peakall and Smouse 2006) with 95% confidence intervals generated from 999 permutations, omitting any populations with fewer than seven samples.

Results and Discussion

AFLP

We found 123 variable loci using the two AFLP primer pairs (Supplementary Data File, AFLP data tab). Of these loci, 106 (86%) were polymorphic at ≥5% level when including all three species. When testing for PCR error in the AFLP process, we found 5 mismatches (i.e., a peak in one run, no peak in the repeat run) in the 48 plants repeated (48 repeats × 123 loci = 5,904 peaks checked for error). This calculates as a 0.08% error rate, which is 0.10 loci in error per plant; thus we considered any samples that were not identical for AFLP genotype as distinct genotypes. Dice pairwise similarities between plants in D. fullonum in the native range varied from 0.56 to 1.00 (identical), and from 0.71 and 0.72 to 1.00 in United States and Argentina, respectively. Native D. laciniatus varied from 0.68 to 1.00, and in the United States varied from 0.83 to 1.00. Dipsacus sativus in the United States (the only collections of this species) varied from 0.90 to 1.00.

Mean population G/N values (only including populations of at least seven individuals) for each of the three species in the United States were significantly different (Likelihood ratio (LR) χ² = 11.37; 2 df; P = 0.003; Figure 3). In a post hoc test, the mean population level G/N value of D. fullonum (0.39) did not differ significantly from D. laciniatus (0.39; odds ratio = 1.09, P = 0.960). However, D. sativus (mean = 0.10) had significantly lower G/N values than both D. fullonum (odds ratio = 4.28, P = 0.010) and D. laciniatus (odds ratio = 3.93, P = 0.030).

Figure 3. Box-and-whisker plot of population G/N values (y axis) for each species (native and invasive samples). D.f., Dipsacus fullonum; D.l., Dipsacus laciniatus; D.s., Dipsacus sativus. Populations with fewer than seven samples not included in analysis.

The proportion of unique genotypes per population also differed in the native versus introduced regions. Populations of D. fullonum in the United States and Argentina (introduced range) did not differ in the mean proportion of unique genotypes (mean G/N = 0.350 [Argentina] vs. 0.324 [United States]; odds ratio = 1.12, P = 0.940). However, both Argentina and the United States had lower G/N values compared with Eurasian populations (native range) of that species (mean G/N [Europe and Asia] = 0.591 vs. Argentina: odds ratio = 0.39, P = 0.013; vs. United States: odds ratio = 2.86, P < 0.0001). Similarly, populations of D. laciniatus in the United States had fewer unique genotypes per population compared with the native Eurasian range (0.395 [United States] vs. 0.734 (Eurasia); LR χ² = 21.70, P < 0.0001) (Figure 4).

Figure 4. Population G/N values for each species and region. Populations with fewer than seven samples not included in the analysis. Letters indicate significant differences between populations in post hoc tests, and the horizontal lines within boxes indicate median values. For the within-species comparisons, the native range (Eurasia) is indicated by a lighter color.

The nuclear ITS region provided 620 bp, of which 26 were variable (4.2%). The haplotype network (Figure 5) contained no homoplasious sites (i.e., no identical mutations found in multiple places on the haplotype network).

Figure 5. Nuclear internal transcribed spacer (ITS) tree for three Dipsacus species; 11 plants sequenced. Hash marks are single-nucleotide changes; boxes indicate genotypes found; and box size indicates relative frequency (D. sativus, n = 3; D. fullonum, n = 6 and n = 1; D. laciniatus, n = 2).

In the STRUCTURE analysis, selection of K for all samples gave a result of K = 2 (Figure 6A). We expected a result of K = 3, given the visual clustering of the PCoA (Figure 7), and suspect that STRUCTURE did not recognize the cluster for D. sativus due to the lower sample size and lower level of variation (most AFLP genotypes were identical or very similar) found in that species. We therefore proceeded with an assumption of K = 3 for the analysis. The STRUCTURE analysis selection of K for D. fullonum native and invasive samples gave a result of K = 2 (Figure 6B).

Figure 6. Delta K result for (A) 572 Dipsacus AFLP genotypes and (B) 361 Dipsacus fullonum AFLP genotypes.

Figure 7. Principal coordinates analysis (PCoA) from Dice similarity data of 572 Dipsacus amplified fragment length polymorphism (AFLP) genotypes. Blue symbols indicate D. fullonum, red indicate D. laciniatus, and green indicate D. sativus.

Most of the genetic variation was found among populations for each species (Table 1), with D. sativus having the highest amount of among-population differentiation due to populations being made up of identical AFLP genotypes. Φ_PT was very similar between U.S./Argentinian/Eurasian D. fullonum (85% to 88%). The species with the highest within-population differentiation was D. laciniatus.

Dice similarity trends (Table 1) show that there is more genetic variation (i.e., plants can have more dissimilar AFLP genotypes) in the native range compared with the invasions, likely indicating a founder effect typically found in invasions or a post-introduction bottleneck (Dlugosch and Parker 2008). Neubert and Caswell (2000) demonstrated that the invasion speed of D. sylvestris (= D. fullonum) was greater than would be expected from demographic models of population increase. Secondary dispersal by different vectors may push range expansion (e.g., Lake et al. 2020) into previously unoccupied areas, and self-compatibility combined with disturbance may strengthen founder effects. The lack of diversity within populations of D. sativus (only one genotype per population) could be attributed to strong founder effects, active selection by humans before naturalization (this is the species historically grown for processing of wool), strong bottlenecks, and/or higher rates of self-pollination than in the other species.

Reproduction

Cross-pollination is noted to be the most common method of reproduction for D. fullonum (Werner 1975) and D. laciniatus (Verlaque 1985), but we found significantly higher G/N values in the native versus invasive range for both species (Figure 4), suggesting higher levels of self-pollination in the invasion compared with their origins. G/N and Φ_PT measurements (Table 1) support that D. sativus has the lowest diversity of the three species and highest amount of among-population differentiation for the three species (Φ_PT in Table 1), and this is likely due to populations being made up of identical AFLP genotypes, suggesting a predominantly self-pollinating reproductive mode. Other possible explanations for D. sativus low invasion diversity are low propagule pressure and resultant inbreeding or our sampling fewer populations of D. sativus than the other two more common species.

Hybridization

In the native range, hybrids have been reported between D. fullonum and D. laciniatus (Gleason and Cronquist 1991; Natural History Museum 2013), but the frequency of these hybrids has not been reported (Gucker 2009). Hybrids are also thought to exist between D. fullonum and D. sativum (Natural History Museum 2013). Werner (1975) notes that plants having intermediate characteristics are found only rarely and that no hybrids have been described or named. We found no heterozygous loci in the nuclear DNA ITS sequences, and thus no indication of recent hybridization. In a review, Solecki (1993) noted D. fullonum and D. laciniatus are only occasionally found together. Our population 58 from Illinois, USA, was morphologically identified as a mix of D. fullonum and D. laciniatus, with all 10 samples (5 of each species) from within a 50-m radius, and we suspected that it would be a highly likely place to find hybrids. The STRUCTURE analysis from that population showed >99% assignment to either species for each plant, and the nuclear DNA had no heterozygous loci, thus there was no indication of hybridization in our collections. We found 6 out of 52 invasive plants with STRUCTURE assignment to a single species at <99% (plant nos. 2, 26, 45, 119, 155, and 402 with percent assignment to species at 87% to 98%; Figure 8 and Supplementary Data File, K = 3 assignment tab), perhaps suggesting some previous gene flow between species, but not recent hybridization (i.e., F₁ hybrids should assign at ∼50% to each paternal species, and backcrosses should assign at ∼75%:25%), although precise assignment of hybrid class can be more complex than stated here (Wringe et al. 2016).

Figure 8. STRUCTURE analysis for K = 3 for 572 plants of Dipsacus fullonum (blue), Dipsacus sativus (green), and Dipsacus laciniatus (red). Bar height within one column (one individual plant) can vary from 0% to 100% assignment value (0 to 1.00 on the y axis). Mixed colors within a column (individual) indicate assignment to multiple species.

Taxonomy

Dipsacus sativus has been named as a subspecies and variety of D. fullonum [Dipsacus fullonum ssp. sativus (L.) Thell. and Dipsacus fullonum var. sativus L.; Missouri Botanical Garden 2023] but is accepted as the separate species D. sativus (L.) Honck. in publications such as Jepson Flora of California (Jepson Flora Project 2023). Our AFLP data showed Dice similarity of ≤48% between D. sativus and D. fullonum, and our ITS DNA sequence data showed 19 single-nucleotide polymorphisms (3.1% sequence divergence) between the two taxa, the same sequence divergence as between D. fullonum and D. laciniatus; thus, both sets of genetic data suggest that D. sativus is a distinct species and not a subspecies or variety of D. fullonum.

Origins

The closest genetic similarities for D. fullonum from United States to the native range were with population 32 in Hungary (Dice pairwise similarity = 0.96) and populations 50 and 51 in Spain at 0.93, and matches to other native samples ranged as low as 0.68. There is support for two genetically distinct clusters of D. fullonum in the STRUCTURE analysis (Figure 6B and indicated by ellipses on Figure 9A). These matches to Hungary and Spain are for U.S. D. fullonum from two different genetic clusters, suggesting two different origins of the U.S. D. fullonum invasion. The closest match for D. fullonum between Argentina and the native range was with Spain population 49 (Dice = 0.94), which is genetically very similar to the same native population that matched with one cluster of the U.S. D. fullonum, suggesting a similar origin from Spain for both the Argentinian and a portion of the U.S. D. fullonum invasions. The next closest country match for Argentina was Greece at 0.87, and values ranged as low as 0.64. Dipsacus laciniatus from the United States most closely matched to population 40 from Russia (Dice = 0.93); the next closest country was Hungary at 0.90, and values ranged as low as 0.72. These highest similarities suggest possible origins of the invasive species. There are cases of host specificity being lower than the species level in biological control programs (Gaskin et al. 2011), and these native locations could be prioritized in searches for potential biological control agents originating from similar plant genotypes. By contrast, an example of a candidate biocontrol agent performing worse on its host population of origin than on different populations of the host plant species has been observed (Cristofaro et al. 2020), highlighting the importance of including multiple populations of a target weed in pre-release evaluations of prospective biocontrol agents.

Figure 9. Principal coordinates analysis (PCoA) from Dice similarity data of (A) 361 Dipsacus fullonum amplified fragment length polymorphism (AFLP) genotypes and (B) 160 D. laciniatus genotypes. Ellipses in (A) indicate U.S. D. fullonum from two different genetic clusters.

Conclusion

In conclusion, levels of diversity and modes of reproduction differ among these three invasive congeners, with D. sativus being lowest in diversity and possibly relying on self-pollination more than the other species; thus it may not have as much potential for evolution of invasive traits or resistance/tolerance to management, though many nondiverse, non-outcrossing terrestrial plant species can be successful, difficult to control invasives (e.g., rush skeletonweed [Chondrilla juncea L.]; Gaskin et al. 2013; Ward et al. 2008). We found no evidence of hybridization within the invasions, though it likely exists outside our collections, and it does not appear to be driving invasion, as occurs in some other species (e.g., Schierenbeck and Ellstrand 2009). Dipsacus sativus is as genetically distinct from D. fullonum as D. fullonum is from D. laciniatus, suggesting that D. sativus is not a subspecies or variety of D. fullonum, and thus may require different management techniques from those applied for D. fullonum. This information regarding invasive teasels’ taxonomy, reproduction, and origins can help us understand their invasive processes as well as give insight into their management.