Hybridization drives evolution of apomicts in Rubus subgenus Rubus: evidence from microsatellite markers

Abstract

Background and AimsRubus subgenus Rubus is a group of mostly apomictic and polyploid species with a complicated taxonomy and history of ongoing hybridization. The only polyploid series with prevailing sexuality is the series Glandulosi, although the apomictic series Discolores and Radula also retain a high degree of sexuality, which is influenced by environmental conditions and/or pollen donors. The aim of this study is to detect sources of genetic variability, determine the origin of apomictic taxa and validate microsatellite markers by cloning and sequencing.

Methods A total of 206 individuals from two central European regions were genotyped for 11 nuclear microsatellite loci and the chloroplast trnL–trnF region. Microsatellite alleles were further sequenced in order to determine the exact repeat number and to detect size homoplasy due to insertions/deletions in flanking regions.

Key Results The results confirm that apomictic microspecies of ser. Radula are derived from crosses between sexual series Glandulosi and apomictic series Discolores, whereby the apomict acts as pollen donor. Each apomictic microspecies is derived from a single distinct genotype differing from the parental taxa, suggesting stabilized clonal reproduction. Intraspecific variation within apomicts is considerably low compared with sexual series Glandulosi, and reflects somatic mutation accumulation. While facultative apomicts produce clonal offspring, sexual species are the conduits of origin for new genetically different apomictic lineages.

Conclusions One of the main driving forces of evolution and speciation in the highly apomictic subgenus Rubus in central Europe is sexuality in the series Glandulosi. Palaeovegetation data suggest that initial hybridizations took place over different time periods in the two studied regions, and that the successful origin and spread of apomictic microspecies of the series Radula took place over several millennia. Additionally, the cloning and sequencing show that standard evaluations of microsatellite repeat numbers underestimate genetic variability considering homoplasy in allele size.

Apomixis, hybridization, microevolution, microsatellites, polyploidy, Rubus subgenus Rubus

INTRODUCTION

The ability to reproduce asexually is relatively common in higher plants, and comprises vegetative reproduction or clonal reproduction through seeds (apomixis or agamospermy), the latter of which is of special interest to geneticists, biosystematists and particularly plant breeders (Barcaccia and Albertini, 2013). Apomixis is found in < 1 % of angiosperms (Whitton et al., 2008) and is mostly missing in crop plants. Although the genetic mechanisms of apomixis have been analysed in several model taxa (Taraxacum, van Dijk and Bakx-Schotman, 2004; Poa, Albertini et al., 2005; Hieracium, Catanach et al., 2006; Hypericum, Schallau et al., 2010), the impact of asexual reproduction on the evolution of many plant taxa remains understudied.

In natural apomictic populations, loss of recombination is hypothesized to lead to (1) deleterious mutation accumulation (Muller, 1964; Lovell et al., 2017; but see Hojsgaard and Hörandl, 2015) and (2) creation of fixed allele combinations resulting in reduced adaptability of clonal genotypes (Kondrashov, 1993), thus leading to the hypothesis that apomixis is an evolutionary dead end. However, highly successful apomicts occur in the families Asteraceae, Rosaceae and Poaceae (Richards, 1986; Whitton et al., 2008). Reasons for their success can be attributed to short-term benefits such as faster colonization (Stebbins, 1950), fixation of successful genotypes (Maynard Smith, 1978) or elevated heterozygosity associated with hybrid origin, polyploidy and asexuality (Lovell et al., 2013). Gene flow between apomicts and sexuals has been documented in many genera (e.g. Pilosella, Krahulcová et al., 2000; Taraxacum, van Dijk, 2003; Majeský et al., 2012, 2015; Rubus, Kollmann et al., 2000; Sochor et al., 2015; Boechera, Lovell et al., 2013; Ranunculus, Hojsgaard et al., 2014) and probably plays an important role in purging asexual genomes from accumulating mutations (Hojsgaard and Hörandl, 2015; Lovell et al., 2017), and producing novel apomictic lineages.

Here we focus on brambles (Rubus subgenus Rubus), whose facultative apomixis and frequent hybridization (Weber, 1996) have led to extensive morphological variability that is reflected in a complex taxonomy consisting of several infrageneric ranks, such as section and series. Out of around 750 species recognized in Europe, sexuals consist of four diploid and two tetraploid (R. caesius, R. ser. Glandulosi) taxa, while all others are polyploids (Kurtto et al., 2010) with varying degrees of pseudogamous aposporous or diplosporous apomixis (Christen, 1950; Pratt and Einset, 1955) with unlinked functional components (i.e. apomeiosis, parthenogenesis and successful endosperm development; Šarhanová et al., 2012). Moreover, three interesting phenomena based on flow cytometric seed screening have been identified: (1) a high degree of residual sexuality in mostly apomictic tetraploid members of the series Discolores and Radula; (2) environmental influences on reproductive mode in R. bifrons (also ser. Discolores); and (3) geographic parthenogenesis in tetraploids of ser. Glandulosi reflected by strict sexuality in the West Carpathians and partial apomixis in the South-west Bohemian Massif (Šarhanová et al., 2012). Whereas ser. Discolores and ser. Glandulosi form two extremes of the Rubus morphological continuum, ser. Radula is intermediate and therefore derived putatively from hybridizations between the two taxa, as is supported by small geographic distribution areas typically characterizing members of ser. Radula, while its putative parents are widespread and thus (expectedly) older (Kurtto et al., 2010).

In the light of previous findings and morphological observations, we focused on members of the series Discolores, Radula and Glandulosi, and asked the following questions. (1) What is the origin of stabilized apomictic microspecies and local hybrids? (2) What is the geographic pattern of inter- and intraspecific variability? (3) To what extent is genetic variation in natural populations influenced by residual sexuality and geographic parthenogenesis? (4) What is the source of inter- and intraspecific variability? Additionally, considering frequent size homoplasy (variability in the sequence of the same size product; Estoup et al., 2002) in microsatellite markers and subsequent errors in population genetic inferences, we validated our microsatellite fragment length analyses by cloning and sequencing products from all loci. The data thus add valuable information on usefulness and caveats of these commonly used molecular markers.

MATERIALS AND METHODS

Plant material

Ten taxa of Rubus subgen. Rubus from two phytogeographic regions of Central Europe were studied: (1) the south-western part of the Bohemian Massif linked geographically to the eastern Alpine region; and (2) the Western (Moravian) Carpathians (n = 206 individuals; Supplementary Data Table S1; Fig. S1). The regions differ by the occurrence of stabilized apomictic microspecies of Rubus ser. Radula, which are missing in the Carpathians (with the exception of R. radula, a widespread and putatively old taxon of different geographic origin which was therefore not included in the study), and by obligate sexuality of the series Glandulosi in the Carpathians compared with facultative apomixis in the Bohemian Massif (Šarhanová et al., 2012). Although many common European brambles are triploid obligate apomicts, we focused on tetraploid members of the subgenus, which maintain the ability to reproduce sexually (Šarhanová et al., 2012; Krahulcová et al., 2013).

Taxa were selected to encompass ongoing evolutionary processes in the studied regions as hypothesized from morphological observations, and included three series from the subgen. Rubus, sect. Rubus. Series Discolores was represented by the tetraploid apomict R. bifrons, which is a common and possibly evolutionarily old taxon with a large distribution range (Kurtto et al., 2010), and by three individuals of the diploid aggregate R. ulmifolius/R. sanctus (hereafter R. ulmifolius agg.), which is a common ancestor to many of the polyploid species of the series and is missing in the studied region (it occurs mainly in the Mediterranean and Atlantic parts of Europe; Sochor et al., 2015). Series Glandulosi was represented by tetraploid accessions without species status (treated as R. hirtus in Kurtto et al., 2010), exhibiting geographical parthenogenesis and prevailing sexuality (Šarhanová et al., 2012). Apomictic ser. Radula was represented by the stabilized microspecies (Dickinson, 1999; Weber, 1996; morphologically uniform biotypes with wide distribution area >50 km, hereafter ‘apo-species’) R. epipsilos, R. indusiatus and R. vatavensis ined. Eight putative F₁R. bifrons × ser. Glandulosi hybrids (as apparent from morphology; Weber, 1996; single bush or single shrubbery, hereafter ‘local hybrids’) producing apomictic seeds were also included. Due to similar hybrid origin (although the sampling was unequal with respect to geography), we merged all of the individuals of these local hybrids into one group. This approach did not affect any result or interpretation of the analyses. Additionally we examined triploid apomict R. grabowskii (one of the most common triploid species in the studied region) and tetraploid apomict R. bohemiicola (very distinct tetraploid morphotype; both from series Discolores). Tetraploid sexual R. caesius from sect. Caesii was selected as outgroup and represents sexual parental species of many members from sect. Corylifolii (Sochor et al., 2015).

Samples were collected from natural populations between 2009 and 2012 (deposited in the herbarium of Palacký University in Olomouc OL), or from the collection of Vojtěch Žíla (deposited in the herbarium of the National Museum in Prague PR and personal collection). To avoid collecting the same ramets (members of one vegetative clone), we collected mostly one individual per species per locality, and considered all localities as populations in major geographic regions.

Chloroplast analyses

DNA was extracted from fresh material or herbarium vouchers following the cetyltrimethylammonium bromide (CTAB) protocol by Doyle and Doyle (1987). The trnL–trnF region was PCR amplified on a GeneAmp PCR System 9700 with primers ‘e’ and ‘f’ (Taberlet et al., 1991). PCR products were purified using the NucleoFast kit (Machery-Nagel), eluted in 30 mL of water and sequenced in both directions on an Applied Biosystems 3730XL sequencer. The sequences were aligned in GENIOUS R6.1.6 (Biomatters). All sequences were deposited in the GenBank nucleotide database (accession nos KX395204–KX395613).

Microsatellite analyses

Forty-three potentially variable microsatellites (Graham et al., 2004; Woodhead et al., 2008, 2013) were chosen for initial tests. PCR of each primer pair was performed on four selected taxa, i.e. R. bifrons, R. epipsilos, R. indusiatus and R. vatavensis, using a QIAGEN Taq DNA Polymerase kit (Supplementary Data Table S2). PCR products were visualized on 3·5 % agarose gels to evaluate potential polymorphisms. For PCR conditions, see Table S2.

Eleven polymorphic loci covering seven linkage groups (Table 1; Woodhead et al., 2008, 2013) were then chosen for fragment size analysis. Forward primers were labelled at the 5′ end with three different fluorescent dyes (Biomers; BMN-5, BMN-6 and DY-751) to allow multiplex analyses. To differentiate clearly peaks belonging to each locus, the PCR was first run separately for each taxon and locus. Allele sizes and intensity of products were determined on the CEQ™ 8000 Genetic Analysis System (Beckmann Coulter) with a 400 bp size standard. To assess reproducibility of the results, random samples were repeated at different stages of analyses.

Table 1.

Descriptive statistics and allelic diversity of all studied microsatellite loci in selected Rubus taxa

Locus	Allele size	No. of repeats	Length of SSR	No. of alleles (size)	No. of alleles (repeats)	LG	Mean size	Mean no. of repeats	Ho	Hs	Private alleles
1B06	192–204	2–5	4	5	4	4	199·1	4	0·918	0·723	–
1G16*	203–247	–	7 + 2	20	–	7	216·5	–	0·971	0·712	Gla-BM, Gla-Ca, ind, vat
72H02	234, 238	10–11	4	2	2	4	234·8	10·2	0·364	0·284	–
Rub238h	140–149	7–10	3 + 3	4	4	3	144·9	8·6	0·862	0·502	Gla × bif, Gla-BM
Rub26a	95–158	18–50	2 + 2	22	21	7	116·2	28·6	1·000	0·795	Gla × bif, Gla-BM, Gla-Ca
Rub35a	246–328	26–67	2 + 2	30	30	5	283·2	44·6	1·000	0·852	Gla × bif, Gla-BM, epi
5K23*	209–231	–	2	11	–	5	223·2	–	1·000	0·707	Gla × bif, Gla-Ca
Rub105b	151–189	7–26	2	16	16	5	169·7	16·4	1·000	0·802	Bif.ca, Gla-BM, Gla-Ca
Rub123a	150–187	12–31	2	16	15	6	160·9	17·9	0·976	0·784	Gla-BM, Gla-Ca, vat
Rub47a	195–245	9–35	2 + 2	18	17	1	213·4	18·3	0·991	0·789	Gla × bif, Gla-BM, Gla-Ca
Rub76b	193–211	2–5	6	4	4	2	201	3·3	0·898	0·634	–
All loci	95–328	2–67		148 (117†)	113				0·907	0·689

No. of repeats, number of observed repeats; no. of alleles, number of alleles counted from the whole amplicon size (size) or number of repeats (repeats); LG, linkage group in the Rubus genome (Woodhead et al., 2008, 2013); H_o, observed heterozygosity; H_s, expected heterozygosity counted from the whole amplicon size;

Taxa and regions abbreviations: Gla, R. ser. Glandulosi; bif, R. bifrons; epi, R. epipsilos; ind, R. indusiatus; vat, R. vatavensis; BM, the Bohemian Massif; Ca, the Carpathians.

*

Loci excluded from analyses based on recognition of the precise number of repeats.

†

Number of alleles without loci excluded from analyses based on recognition of the precise number of repeats.

Visualization and initial determination of correct peak size was done using the CEQ™ 8000 Genetic Analysis Software (Beckmann Coulter) and manually verified by eye to avoid errors due to band stuttering. The observed sizes were compared with results of cloning and sequencing (see below), and sizes of loci with ≥3 bp simple sequence repeat (SSR) motifs were converted to size (bp) from sequencing results. Due to slippage during PCR, it was not possible to determine amplicon size unambiguously in loci with 2 bp motifs, and thus fragment size was corrected for the number of repeats divided by two.

Cloning and sequencing of microsatellite amplicons

For all SSR loci in each taxon, PCR amplicons were cloned using a CloneJET PCR Cloning Kit (Fermentas) following the manufacturer’s protocol, and transformed into DH5a Escherichia coli strains. Positive colonies were transferred to 200 mL of LB liquid medium with 0·1 mg mL^–1 ampicillin and incubated overnight at 37 ºC, and 1 μL of product was used for unidirectional sequencing using pJET-F and pJET-R primers. The number of sequenced colonies per locus varied according to the number of alleles (based on amplicon size) per sample – four, if only one allele was detected, eight if 2–5 alleles were present. The procedure was repeated for selected loci to detect more variants (Supplementary Data Table S3). If not all sizes (as seen in fragment analyses) were recovered during sequencing, PCR followed by cutting specific fragments from an agarose gel, purification with a QIAquick PCR Purification Kit (Qiagen) and cloning of the selected fragments was performed. All sequences of SSRs and their flanking regions (FRs) were deposited in the GenBank nucleotide database (accession nos KX395204–KX395613).

As allele cloning and sequencing were performed to determine the origin of size variability in SSRs and FRs, any single nucleotide polymorphisms (SNPs) detected in FRs were not considered in subsequent analyses. To calculate the real number of repeats for each locus we used the formula:

[whole amplicon size (bp) – the smallest observed size of FRs – primer sizes]/length of repeat motif

Statistical analyses of SSR data

Analyses were based on (1) total length of the fragments at all 11 loci and (2) number of SSR repeats based on sequencing of nine selected loci (Loci 1G16 and 5K23 did not fit the stepwise mutation model due to variation in repeat motifs and therefore were excluded).

SPAGEDI v1.4c (Hardy and Vekemans, 2002) was used to calculate the total number of alleles/repeats, mean allele size/number of repeats and range of allele sizes/number of repeats. Global F- and R-statistics were calculated for both data sets separately (see above), in addition to a data set created by dividing the total amplicon size by the length of the SSR motif. To estimate the effects of the infinite allele (IAM; Kimura and Crow, 1964) vs. stepwise mutations (SMM; Ohta and Kimura, 1973) models, we compared F_ST vs. R_ST and tested whether observed R_ST is significantly larger than its value after permutation. P-values were obtained after 10 000 permutations of allele sizes among alleles within loci to test the null hypothesis that ‘stepwise mutations do not contribute to genetic differentiation’ (Hardy et al., 2003).

Genotypic variation was determined in GENODIVE 2.0b25 (Meirmans and van Tienderen, 2004). We calculated multilocus genetic distances from the number of repeats at nine loci using the SMM model (based on recognition of the precise number of repeats) and from the fragment size of all loci using the IAM model. The histograms of frequency distribution of all pairwise distances were used for threshold determination to distinguish between inter- and intraclonal diversity, whereby closely related individuals (below a threshold) are considered to be the same clone. The histograms were also generated for each taxon separately using all loci and IAM. GENODIVE identified genotypic diversity and evenness, and F_ST values for all pairs of populations. Mutation and recombination as sources of intraspecific variation were also distinguished by the character incompatibility approach in module JACTAX of package PICA 4.0 (Wilkinson, 2001), which is based on the rationale that binary unordered characters (i.e. presence or absence of alleles) cannot be present in all combinations of character states in a population in the absence of recombination. If the characters are compatible (i.e. two binary characters form only 1–3 combinations), mutations are a very probable source of the observed variation (Mes, 1998).

Pairwise genetic identities between all individuals were calculated from multilocus data of the total fragment size in POPDIST (Guldbrandtsen et al., 2000) with the Tomiuk and Loeschcke (1991) algorithm, which is suitable for comparison of genotypes with different ploidy levels and recommended for asexual taxa. The matrix was used to generate a principal co-ordinates analysis (PCoA) in SPSS v20 (release 2011, IBM Corp.). Additionally we performed a Mantel test to detect association between genetic similarity and geographic distances of each taxa and each region using the Isolation by Distance web service v3.23 (Jensen et al., 2005) with 10 000 randomizations.

RESULTS

Series Radula and Glandulosi share the same trnL–trnF haplotype

The chloroplast trnL–trnF region of 193 individuals showed six haplotypes (Table S1), and was 450 bp long for R. caesius and 444 bp for all other taxa, with the latter haplotypes differing by single substitutions (Fig. 1; Supplementary Data Table S4). All individuals of R. bifrons shared a single species-specific haplotype.

Fig. 1.

Haplotype network derived from sequences of the chloroplast trnL–trnF region of different taxa of the subgen. Rubus, using R. caesius as an outgroup. Black dots indicate the number of mutational changes, species belonging to ser. Radula are underlined, and species belonging to ser. Discolores are in bold.

Sequencing of microsatellites allows unambiguous allele assessment and reveals common homoplasy

Sequencing of SSR amplicons (Supplementary Data Table S3; Fig. S2) enabled the elucidation of ambiguous fragment sizes, as were identified for loci 5K23 (168, 176, 246 and 262 bp), Rub26a (94, 96 and >140 bp) and Rub35a (>300 bp; for full names of the loci, see Table S2). Peak heights of Rub76b were highly variable within individuals, with some peaks just above the detection level (for examples of capillary sequencing electropherograms; see Supplementary Data Fig. S3). All ambiguous products of 5K23 were not detected, even after sequencing of isolated bands, and were hence discarded. In contrast, ambiguous alleles of loci Rub26a and Rub35a, in addition to small peaks of Rub76b, were all confirmed by sequencing.

Sequencing enabled us to calculate precisely the number of repeats for most loci. Locus 72H02 contained repeat variability, which differed from its published sequence (Woodhead et al., 2008; TAAA vs. CTAG, the latter of which was also detected within the SSR sequence but lacked any variability), while loci Rub238h and 1G16 had additional variable repeat regions (ACG and TTAAAAA, respectively). Loci with 2 bp motifs had highly variable sizes as recognized via sequencing, which were also observed from fragment analyses as stutter peaks. At five loci (Rub26a, Rub35a, Rub47a, 1G16 and Rub238h), two independently derived SSR motifs were detected. For locus 1G16, the number of repeats was not possible to determine from the amplicon size as it was characterized by two variable repeat motifs of different size (7 and 2 bp). Locus 5K23 was highly variable due to insertions within the repetitive region. Both loci were excluded from analyses of repeat numbers.

We detected a high number of SNPs in most repeats and FRs which mostly did not correspond to variation in the number of SSR repeats, which led to size homoplasy at most loci (Fig. S2).

In 5K23 and Rub123a, a single base pair indel in the FR led to different alleles, and in locus 1G16 a single base pair difference was caused by two different repetitive motifs, 7 and 2 bp in size. Thus these double peaks observed in the fragment analyses were accepted as different alleles for all individuals (see an example of the electropherograms in Fig. S3). In cases where we characterized alleles at these loci by the number of motif repeats, the 1 bp difference of the flanking region was lost. Hence, as we were interested in examining whether the 1 bp variation enabled greater resolving power in our attempt to identify different apomictic lineages, we created two data sets based upon (1) the presence of each allele as defined by its length including the 1 bp variation and (2) the allelic state characterized by the number of repeats only (i.e. without considering the 1 bp variation).

Genotypic variation reflects reproduction mode

Out of 11 microsatellite loci, five were highly polymorphic, with >20 alleles each (Table 1; Supplementary Data Table S5). In contrast, four microsatellites characterized by repeat motifs >2 bp had few polymorphisms and less than six alleles, suggesting a negative correlation between polymorphism and repeat motif size. In total, 148 alleles were identified among loci based on amplicon size, and 113 alleles from nine selected loci based on the number of repeats (see the Materials and Methods; Table 1). Due to ignoring 1 bp indels in calculations of the number of repeats the number of alleles dropped from 117 to 113 among the nine selected loci compared with the total amplicon size. This resulted in a considerable decrease in the number of all genotypes, especially for apomictic taxa (Table 2). Several tetraploid individuals had more alleles than expected as compared with their ploidy, and this was evident in eight individuals of R. epipsilos at locus Rub123a, two individuals of ser. Glandulosi at 5K23 and two other individuals of ser. Glandulosi at 1B06. In the statistical analyses all these alleles were treated as present. Private alleles belonged mainly to sexual ser. Glandulosi and hybrid R. bifrons × ser. Glandulosi, although several private alleles (1–3) occurred in all other apomictic taxa (R. bifrons, R. epipsilos, R. indusiatus and R. vatavensis; Tables 1 and 3).

Table 2.

Summary of multilocus descriptive statistics of selected tetraploid populations of subgen. Rubus (taxa names abbreviated as in Table 1) based on 11 loci counted from the total fragment size and IAM model (size), and nine loci counted from the exact number of repeats and SMM model (repeats)

	bif-BM	bif-Ca	epi	ind	vat	Gla × bif	Gla-BM	Gla-Ca
Sample size	14	29	36	46	26	8	19	14
tNG, size/repeats	11/6	14/6	26/15	14/3	14/6	8/8	19/18	14/14
NG, size/repeats	1/1	1/1	1/1	1/1	1/2	7/5	15/11	14/14
eve, size/repeats	1·0/1·0	1·0/1·0	1·0/1·0	1·0/1·0	1·0/0·54	0·914/0·64	0·832/0·67	1·0/1·0
Div, size/repeats	0·0/0·0	0·0/0·0	0·0/0·0	0·0/0·0	0·0/0·077	0·964/0·786	0·965/0·912	1·0/1·0
Ho, size/repeats	0·957/0·915	0·96/0·918	0·953/0·906	0·909/0·852	0·944/0·915	0·951/0·88	0·809/0·752	0·776/0·725
Hs. size/repeats	0·663/0·65	0·653/0·641	0·667/0·658	0·624/0·579	0·651/0·647	0·784/0·741	0·738/0·685	0·74/0·705
GIS, size/repeats	–0·444/–0·409	–0·471/–0·434	–0·43/–0·378	–0·458/–0·472	–0·449/–0·414	–0·213/–0·188	–0·096/–0·098	–0·049/–0·028
MIC, size/repeats	8/4	6/7	19/19	6/6	17/15	563/347	1404/791	1432/946

apo, apomictic taxa; sex, sexual taxa; tNG, total number of all genotypes; NG, number of recombinant genotypes; eve, evenness (= effective number of genotypes/tNG); div, Nei’s (1978) genotypic diversity; H_o, observed heterozygosity; H_s, expected heterozygosity; G_IS, inbreeding coefficient; MIC, matrix incompatibility count.

Table 3.

Shared alleles among the studied taxa

	R. bifrons		R. ser. Radula			R. ser. Glandulosi
	bif-Ca	bif-BM	epi	ind	vat	Gla-BM	Gla-Ca
Rub47a	201/211/217/219	201/211/217/219	201/211/223/245	201/211/230	195/201/211	195/201/211/217/219/223/230/243	195/201/211/217/219/223/243/245
Rub76b	199/205	199/205	193/199/205	199/205/211	193/199/205	193/199/205/211	193/199/205
Rub105b	171/173/175/179	171/173/177/179	157/159/171/179	161/165/171/177	167/173/177	157/159/161/165/171/173/175/177	157/159/161/165/167/171/173/175/179
Rub123a	152/153/167/175	152/153/167/175	152/153/157/159	153/157/167	153/167/183 (185)	153/157/159/167/175	152/153/157/159/167/175/185
5K23	223/228	223/228	221/223/228	209/222/223/228	221/223/228	209/221/222/223/228	209/221/222/223
Rub238h	140/149	140/149	140/149	140/149	140/149	140/149	140/149
72H02	234/238	234/238	234/238	234	234/238	234	234/238
1G16	203/208/216	203/208/216	208/221	203/216/230	216/245	203/216/221/232/237	203/216/221/233/236
1B06	198/200/204	198/200/204	198/200/204	192/198/200	196/198/200/204	192/196/198/200/204	192/196/198/200/204
Rub26a	113/115/117	113/115/117	95/115/117/119	111/113/117/141	111/113/115/117	97/113/115/117/119/141	97/111/113/115/117/119/141
Rub35a	268/272/318/322	268/272/318/322	246/272/284/322	268/272/316	260/268/272/322	268/272/284/318	268/272/284/318

Alleles not found in ser. Glandulosi are underlined, not found in R. bifrons are in bold, and alleles specific for ser. Radula are in italics. Primary hybrids and private alleles of R. bifrons and ser. Glandulosi are not shown for clarity. For the full data set, see Supplementary Data Table S9. Taxa names are abbreviated as in Table 1.

To elucidate the origin of genetic variability in each taxon, we used two different approaches. First, we created histograms of pairwise genetic distances in GENODIVE, with the smallest frequency class serving as a threshold which differentiated between variability arising via sexual reproduction vs. mutation accumulation within an apomictic clonal lineage (Supplementary Data Fig. S4). Secondly, the number of marker incompatibility conflicts was counted in JACTAX. When the total amplicon length and IAM were applied, we detected recombinant genotypes (genetic variability resulting from sexual reproduction) only for ser. Glandulosi [mean genetic distance 13, matrix incompatibility count (MIC) = 1404 for the Bohemian Massif and 1432 for the Carpathians] and local hybrid individuals R. bifrons × ser. Glandulosi (mean genetic distance 13, MIC = 563). Within the former, three genotypes were identified which were represented by three, two and two individuals, respectively, all from the Bohemian Massif (Supplementary Data Table S6). When the precise number of repeats and SMM were used, an additional genotype was observed for R. vatavensis, though the diversity was very low (0·077; Table 2) whereby the second genotype was a single individual characterized by a private allele at locus Rub123a which caused only four conflicts in the incompatibility analysis. In contrast, histograms of genotypic distances with a single peak around zero indicated no evidence of sexually derived variation in other taxa (Fig. S4). The MIC was also low for each of the taxa (not exceeding 19; Table 2), as well as the contributions of each individual (not exceeding eight conflicts; data not shown).

Pairwise F_ST and R_ST indices showed that genetic distances between every two species/regions were roughly similar (Supplementary Data Table S7). The R_ST value was not significantly larger than R_ST after permutation, which indicated no contribution of pairwise mutations to genetic differentiation (Supplementary Data Table S8), and thus the IAM model was considered for all studied loci. Individual inbreeding coefficients (G_IS;Table 1) varied from –0·028 to –0·098 for sexual ser. Glandulosi, and from –0·378 to –0·472 for the apomicts, reflecting clonal diversity. The value was intermediate for local hybrids R. bifrons × ser. Glandulosi (–0·213 or –0·188 based on total amplicon size or number of repeats, respectively). Observed multilocus heterozygosity (H_o) was higher than the expected heterozygosity (H_s) within populations of apo-species and also in local hybrids, whereas both values were similar for sexuals (Table 2).

The PCoA (Fig. 2) clearly differentiated all morphologically distinct apo-species. The distribution of sexually reproducing biotypes of ser. Glandulosi reflected high genotypic variability, and both the local hybrids R. bifrons × ser. Glandulosi and apo-species from ser. Radula appeared intermediate on the first axis between ser. Glandulosi and the tight cluster of R. bifrons. A mantel test (Supplementary Data Fig. S5) for correlation between genetic and geographic similarities showed significant isolation by distance (IBD) for R. bifrons (r = –0·582, P < 0·0001) and ser. Glandulosi (r = –0·1033, P = 0·0154) when analysed for the whole study area. The correlations for each region separately remained with low significance for ser. Glandulosi in the Carpathians (r = –0·1787, P = 0·0461) and additionally for R. vatavensis occurring only in the Bohemian Massif (r = –0·4200, P = 0·0065).

Fig. 2.

Principal co-ordinate analysis based on Tomiuk–Loeschcke (1991) identities of 11 loci of all studied Rubus taxa. The first two axes explained 72·6 % variability. Ellipses delimit taxa involved in hybridization processes. bm, Bohemian Massif; ca, the Carpathians.

DISCUSSION

Here we have used microsatellite fragment analyses of 11 SSR loci followed by cloning and sequencing to validate the data for SSR repeat polymorphism analysis. Together with chloroplast trnL–trnF sequences, we were able to elucidate the origin of several apomictic species of European brambles and assess the evolutionary mechanisms contributing to success of apomictic taxa of R. subgen. Rubus in Central Europe.

Homoplasy is a common phenomenon

As with neutral and repetitive sequences, microsatellites exhibit elevated mutation rates which increase the probability that two alleles identical by state (sequence and/or length) are not identical by descent (i.e. homoplasy; Estoup et al., 2002). Nucleotide substitution in these regions is believed to evolve more slowly than the number of repeats of the microsatellite itself, but the flanking regions can exhibit some degree of polymorphism, and in some cases significantly (Väli et al., 2008; Corral et al., 2009). Since the amplified PCR fragments may include diverse mutations in FRs and/or SSRs which sometimes contain multiple repetitive motifs (Fig. S2), homoplasy may be a serious problem when only fragment length analysis is performed, leading to underestimation of genetic variability and consequently to errors in population genetic inferences. Although the chance of homoplasy is proportional to the genetic distance between two individuals (Selkoe and Toonen, 2006), the few existing sequence-based studies revealed a high degree of fragment length homoplasy not only among different species but also within one species or a population (Estoup et al., 2002; Curtu et al., 2004; Anmarkrud et al., 2008; Javed et al., 2014; Germain-Aubrey et al., 2016). Our sequence analysis provide clear evidence of homoplasy at loci containing (1) SNPs or indels in FRs (5K23, Rub238h, Rub123a, 1G16, 1B06, 72H02 and Rub26a); (2) multiple repeat motifs (compound SSR; Rub47a, Rub238h, 1G16, Rub26a and Rub35a); and (3) interrupted (imperfect) repeats (mainly 5K23). Moreover, if we consider SNPs as evolutionarily more stable and more informative than SSRs (Estoup et al., 2002), discrepancies between SNPs and SSRs at several loci could also be explained by convergent evolution in allele length (see Fig. S2). Nevertheless, our assessment of clonal diversity was not markedly affected by homoplasy considering high allelic variation and thus good discriminatory power of microsatellites (Estoup et al., 2002).

In most studies the number of SSR repeats is calculated by dividing the whole amplicon size by the repeat motif length, a step which can be misleading since invariable FRs are assumed. However, as shown in our study (Supplementary Data Fig. S2; Table S3) and known from others (Mogg et al., 2002; Curtu et al., 2004; Germain-Aubrey et al., 2016), this is often not the case, and is especially prone to error in apomicts which, by definition, have disturbed meiosis leading to accumulation of mutationa (Corral et al., 2009). Sequencing of SSR loci can help to determine the real number of repeats and/or additionally confounding SNPs and indels. For example, 1 bp variability is overlooked in most studies, although it can be easily detected by sequencing (see examples of electropherograms in Fig. S3). We detected such single base pair variability in three loci out of 11, suggesting that it is not a rare phenomenon. Disregarding such information leads to an underestimation of the number of genotypes (tNG and NG; Table 1), especially for apomictic taxa.

Intraspecific variability arises via mutation, and confirms the species concept

The three main sources of genetic variability of apomictic plants are (1) the origin of the apomicts from different sexual backgrounds, (2) mutation accumulation (i.e. Muller’s ratchet; Muller, 1964) and (3) hybridizaton of apomictic and sexual genotypes, which together yield two different patterns of variability. Of these, variability driven by mutations is lowest and by definition has had insufficient time to accumulate in relatively young apomictic lineages (or taxa). In general, facultative aposporous apomicts retain the ability to reproduce sexually and to produce viable pollen (Bicknell and Koltunow, 2004), and hence their genotypic variability is often comparable with their sexual relatives (e.g. Ranunculus kuepferi, Cosendai et al., 2011, 2013; Ranunculus variabilis, Hörandl et al., 2001; Erigeron compositus, Noyes and Soltis, 1996; reviewed by Hörandl and Paun, 2007). This is also the case in the subgen. Rubus where both processes (sexual and apomictic) can take a place within one individual or fruit. Normal sexual development (i.e. meiotic reduction of the megaspore and subsequent fertilization) was observed in all tetraploid apomictic species studied here, with frequencies ranging from 9 % in Bohemian R. bifrons to 23 % in R. epipsilos, the proportion of which may even increase under specific environmental conditions (Šarhanová et al., 2012). Moreover tetraploid species mostly exhibit high pollen viability (Nybom, 1988), further increasing the probability of recombination.

In contrast to the high degree of sexuality in apomictic species, we detected no traces of recombination on the intraspecific level in R. epipsilos and R. indusiatus, and thus the observed variation has probaby been generated by mutations (Table 2). The slightly increased variability of R. vatavensis (compared with R. epipsilos and R. indusiatus) can also be attributed to mutations, although multiple origins from genetically similar individuals or ongoing recombination cannot be excluded. High interspecific variation reflects an independent origin of each apo-species from a different sexual background (Fig. 2), a pattern ultimately associated with the species concept (Weber, 1996; Holub, 1997) and precise clone identification based on morphological observations. This is in contrast to many other apomictic taxa with multiple origins (e.g. Crateagus, Talent and Dickinson, 2005; Paspalum, Siena et al., 2008; Taraxacum, van der Hulst et al., 2000; Ranunculus kuepferi, Cosendai et al., 2011). Similar to this study, Majeský et al. (2012, 2015) demonstrated low genetic variability and single hybrid origins of different apomictic microspecies in central European Taraxacum. Nonetheless, discussion on the microspecies concept in Rubus is beyond the scope of this work and here we simply accept the current view (Weber, 1996; Holub, 1997) by confirmation of (formerly only supposed) intraspecific genetic homogeneity.

Residual sexuality drives bramble evolution

In most apomictic complexes, the sympatric occurrence of sexuals and apomicts facilitates hybridization, although they often differ in ploidy levels – the sexuals being diploid and asexuals being polyploid (Asker and Jerling, 1992). Thus hybridization can sometimes be blocked by incompatibility of ploidy levels or endosperm imbalance of maternal vs. paternal genomes (Nogler, 1984; Vinkenoog and Scott, 2001). The latter does not seem limiting for brambles, as they are characterized by relaxed pressure to maintain the normal 2m:1p (maternal genome:paternal genome) endosperm balance number, which can vary up to 8m:1p or occasionally reflect autonomous development (i.e. 8m:0p; Šarhanová et al., 2012). This is possible considering that apomictic brambles also retain the ability to develop viable pollen of different ploidies (Nybom, 1988), which can fertilize both sexual and facultatively apomictic individuals.

Despite limited possibilities to visualize relationships of allopolyploid species and the high number of shared alleles, our data confirm the morphological observation that stabilized apo-species and local hybrids of ser. Radula were formed through a cross between tetraploid sexual ser. Glandulosi and tetraploid apomictic R. bifrons, or other closely related (extinct or unsampled) apomictic taxa (Table 3; Fig. 2). The hybrid origin is also supported by genome size variation, where individuals of ser. Glandulosi have smaller genomes compared with R. bifrons, while the studied species from ser. Radula have an intermediate genome size (Šarhanová et al., 2012).

The co-occurrence of polyploid sexuals and apomicts has led to the creation of hybrid swarms in Hieracium subgen. Pilosella (Krahulcová et al., 2009), where apomicts can serve as pollen donors and seed parents (Fehrer et al., 2007; Krahulec et al., 2008) and the residual sexuality of apomicts leads to the creation of new variability and formation of new apomictic lineages (Houliston and Chapman, 2004; Morgan-Richards et al., 2004). These phenomena were also observed for R. subgen. Rubus, particularly in the sect. Corylifolii, which contains numerous apo-species and local biotypes which are all derived from hybridizations between tetraploid sexual R. caesius (serving both as a seed and a pollen parent) and other brambles (Sochor et al., 2015).

In contrast to the above-mentioned studies on hawkweeds and different brambles (Fehrer et al., 2007; Krahulec et al., 2008; Sochor et al., 2015), our chloroplast data point to unidirectional hybridization between apomicts as pollen donors and sexuals as mother plants (Fig. 1; Table S4). A selective process can be hypothesized which prevents hybridization with an apomictic maternal background, although it is unclear whether this involves a pre- or post-zygotic mechanism. Alternatively, lower proportions of reduced megaspores in apomicts (Šarhanová et al., 2012) may simply decrease the probability of an apomict being a seed parent in the cross, and thus unidirectional hybridization then would not be an absolute rule. This explanation supports the study by Sochor et al. (2015) of R. radula (a widespread species from ser. Radula) which was derived from the reciprocal cross ser. Discolores × ser. Glandulosi.

Geographic parthenogenesis may be shaped by selection

Geographic parthenogenesis refers to the occurrence of apomicts in larger (or different) areas compared with their sexual relatives (Vandel, 1928), and is known for many agamic complexes (e.g. Ranunculus, Hörandl, 2009; Cosendai et al., 2013; Taraxacum, van Dijk, 2003; Crataegus, Lo et al., 2012), although some authors pointed out that the difference in geography may be primarily caused by different ploidy levels (Bierzychudek, 1987; van Dijk, 2003). Similarly, ecological niche shift was explained by polyploidization, rather than reproduction mode in Boechera (Mau et al., 2015). Geographic parthenogenesis is well known in European brambles on the subgenus level (Asker and Jerling, 1992), and was further detected by flow cytometry seed screening of the tetraploid ser. Glandulosi which is strictly sexual in the Western Carpathians and partly apomictic in the South-west Bohemian Massif (Šarhanová et al., 2012). Here, no evidence of clonal genotypes was found in the Carpathians, while three genotypes represented by three, two and two individuals, respectively, were detected in the Bohemian Massif. Accordingly, genotypic diversity was slightly lower and observed heterozygosity slightly higher in the latter region (Table 2). Although we detected IBD (Supplementary Data Fig. S5) for both ser. Glandulosi and R. bifrons when considering the whole studied area, PCoA (Fig. 2) or any other clustering method did not show any pattern of genetic structuring among the regions based on the neutral SSR markers, nor did they detect any differentiation between facultatively apomictic and sexual individuals of ser. Glandulosi (data not shown). IBD in ser. Glandulosi may thus reflect reduced gene flow between two geographically isolated metapopulations, rather than differences in reproduction mode. In the case of R. bifrons, the observed divergence was caused by a single allele mutation (175 vs. 177 on the Rub105b locus; Table 3), indicating differentiation by mutation among the two lineages rather than recent gene flow between them.

Without further study, causes of geographic parthenogenesis remain unclear. Apomictic organisms may be more successful than their sexual counterparts at margins of occurrence and through their ability to spread from single individuals with a fixed successful genotype (Maynard Smith, 1978; Hörandl, 2006). This advantage may be reduced in pseudogamous apomicts, if fertilization of the central cell via outcrossing is necessary for endosperm development, although this is no handicap for subgen. Rubus in which endosperm can develop via selfing (Nybom, 1988). Nevertheless, Haskell (1960) observed a mentor effect whereby apomictic brambles have better fitness and higher germination frequencies when endosperm originates via outcrossing rather than selfing, and increased germination when ‘outcrossed’ with higher ploidy pollen donors. Kollmann et al. (2000) also documented higher seed set after ‘outcrossing’ in apomictic brambles. Assuming similar increases in fitness in natural populations, the availability of outcrossed pollen in brambles is not limiting, as numerous Rubus species and genotypes often co-occur, and the two studied regions do not differ in availability of suitable pollen donors (unpubl. data). The ability to spread from a single individual therefore does not seem to be affected by pseudogamy, and cannot be rejected as a factor potentially shaping geographic parthenogenesis.

Apomictic taxa are expected to have increased heterozygosity compared with their sexual relatives due to their hybrid origin, polyploidy and fixed allele composition through apomixis, which can further increase by mutation accumulation through time (Birky, 1996; Hörandl and Paun, 2007). Alternatively, elevated heterozygosity in apomicts relative to their sexual progenitors could reflect a selective process whereby beneficial heterosis-like effects are maintained in a fixed asexual genotype. Such a pattern was observed for different apomicts (e.g. Boechera, Beck et al., 2012; Lovell et al., 2013; Ranunculus, Paun et al., 2006; Taraxacum, Hughes and Richards, 1988) and was confirmed here (Table 2). Hence, apomixis could be favoured in areas where heterozygosity of sexual populations tends to decrease due to genetic drift or inbreeding, e.g. in small isolated populations (or) at margins of the species distribution area. This hypothesis is supported by our data as the Bohemian Massif population of ser. Glandulosi exhibits higher heterozygosity (Table 2) despite its position at the distribution limit (Kurtto et al., 2010).

Bramble speciation is probably driven by human

Although species diversity of hybridogenous R. ser. Radula in the Southern Bohemian Massif is striking (Kurtto et al., 2010), it is not the only region where both its parental taxa (ser. Glandulosi and R. bifrons) are found in sympatry with their local hybrids. In other areas, the hybrids can be formed just as frequently, for example in the Carpathians where stabilized apo-species of ser. Radula are practically missing (Dančák et al., 2005; Kurtto et al., 2010). Explanations for this pattern may be manifold; for example, we previously speculated that introgression from apomictic R. bifrons may decrease the sexuality of ser. Glandulosi in the Bohemian Massif (Šarhanová et al., 2012). This hypothesis is partly supported by the SSR data since all but two alleles of R. bifrons were found in ser. Glandulosi (Table 3). It is nevertheless contradicted (1) by the absence of marked geographic genetic structure, (2) by three shared SSR alleles between both taxa detected only in the Carpathians and not the Bohemian Massif in ser. Glandulosi (Table 3) and (3) by apomixis in the North Bohemian ser. Glandulosi with an absence of R. bifrons and related taxa (unpubl. data). Also, Sochor et al. (2015) detected high proportions of Glandulosi-like internal transcribed spacer (ITS) ribotypes in the R. bifrons genome, which points rather to shared ancestry and incomplete lineage sorting resulting in high numbers of common SSR alleles. The reverse association between geographic parthenogenesis and apo-species diversity, i.e. that the spread of hybrids is caused by increased apomixis in ser. Glandulosi, also seems unlikely as the hybrids are able to produce apomictic seeds in both studied regions (Šarhanová et al., 2012). Stabilization of the hybrid state may therefore be related to external factors, such as environmental conditions and/or history.

Rubus evolution in Europe has been affected by Holocene migration from glacial refugia (Sochor et al., 2015) and human impact on landscape, with only a limited number of species surviving under (semi-)closed forest canopy. Both hybridization events and apomictic spread of hybrids have thus been facilitated by human-driven deforestation since the Neolithic Era (Matzke-Hajek, 1997). Whereas brambles of ser. Glandulosi are typical representatives of spruce or beech–fir mountain forests, R. bifrons is a light-demanding species of forest edges and clearings, which avoids both forest and open natural habitats and exhibits strong association with human activities (Weber, 1995). Both taxa probably survived the last glacial maximum in the Northern Balkans or the Eastern Alpine Foreland (or evolved there from sexual ancestors) and migrated northwards during the Holocene (Sochor et al., 2015; see also Kurtto et al., 2010). Migration of ser. Glandulosi into both studied regions is expected to have occurred during expansion of spruce, beech and/or fir in the Boreal/Atlantic period between approx. 9·0 and 6·3 cal. ka BP (Svobodová et al., 2001; Jankovská, 2006; Petr et al., 2013). Rubus bifrons could have been present in the Southern Bohemian Massif as early as the later Atlantic (6 cal. ka BP) and even in the mountainous areas (Svobodová et al., 2001; Jankovská, 2006), as palynological evidence for extensive pastures and agriculture suggest. On the other hand, it is unlikely (based on vegetation reconstructions) that this species had appeared in (sub-)mountain areas of the Western Carpathians before the first modern human settlements in the 13th century AD. Its spread in the Carpathians may have occurred even later, possibly in connection with so-called Walachian colonization 600 years ago (Rybníček and Rybníčková, 2008). Thus, the first hybridization events in the Bohemian Massif may have pre-dated those in the Carpathians by even more than five millennia. This time span could have been crucial for the successful origin and spread of apomictic hybrids.

Conclusion

It has been previously observed that the pollen donor is crucial for switching of reproductive mode between apomictic and sexual Rubus (Šarhanová et al., 2012), where absence of outcrossing leads to increased sexuality in apomictic species. Our molecular data show that apomictic brambles are highly heterozygous and that residual sexuality of apomicts does not give origin to successful offspring, and that sexual pollination of apomictic maternal Rubus species is unlikely. Unidirectional hybridization from apomicts to sexuals may lead to the origin of new local (apomictic) hybrids with new genotypes, and may provide a mechanism to escape Muller’s ratchet (well discussed for Taraxacum and Chondrilla; van Dijk, 2003). Together with the results presented here, we hypothesize that facultative sexuality has a negligible influence on the origin of variability of apomictic brambles, but rather the driving force of evolution and speciation of subgen. Rubus in central Europe is the sexuality of ser. Glandulosi (and the other sexual tetraploid, R. caesius; Sochor et al., 2015). It seems that tetraploid apomictic brambles are quite successful and old entities in Europe and have escaped the expected dead end of deleterious mutation accumulation (Maynard Smith, 1978; Kondrashov, 1994) via their ability to create new biotypes through hybridization. Holocene anthropogenic changes in landscape probably enabled both migration of parental taxa and spread of apomictic hybrids, and represents a leading factor of bramble evolution in Europe.

SUPPLEMENTARY DATA

Supplementary data are available online at https://academic.oup.com/aob and consist of the following. Figure S1: sampling design of studied tetraploid taxa of Rubus subgen. Rubus. Figure S2: alignments of SSR sequences sorted by length. Figure S3: examples of fragment analysis of microsatellite loci of different taxa of subgen. Rubus. Figure S4: histograms of the frequency distribution of the genetic distances based on amplicon size of 11 loci and IAM. Figure S5: relationship between Tomiuk–Loeschcke identities (genetic similarity) and geographic distances among individuals of different Rubus taxa using Mantel test with a reduced major axis regression line and probabilities of correlation values. Table S1: geographic co-ordinates and sampling locations of all studied individuals of subgen. Rubus with region of origin, taxa abbreviation and chloroplast haplotype specification. Table S2: conditions of all PCRs used for genetic characterization of selected taxa of subgen. Rubus. Table S3: sequencing results of 11 microsatellite loci of selected taxa subgen. Rubus. Table S4: variable positions in sequences of the trnL–trnF region of chloroplasts distinguishing six haplotypes of all studied species of subgen. Rubus. Table S5: descriptive single locus statistics of eight selected tetraploid populations of subgen. Rubus based on 11 microsatellite loci. Table S6: assigned genotypes of each individual based on absolute allele size and thresholds of 0, 5 and 6 mutations. Table S7: single locus estimates of global F- and R-statistics for selected Rubus taxa counted from whole amplicon size, number of repeats and precise number of repeats estimation via sequencing. Table S8: values of F_ST genetic differentiation between selected tetraploid populations of subgen. Rubus based on total fragment size of 11 microsatellite loci and number of repeats of nine loci. Table S9: observed alleles of all individuals of Rubus subgen. Rubus.

ACKNOWLEDGEMENTS

We thank all colleagues who helped with sample collection, namely M. Lepší, P. Lepší, V. Žíla and A. Jírová. We thank T. Pěnkavová for help with cultivation and P. Oswald for help with cloning experiments, and thank in particular F. Blattner for careful reading of and comments on the manuscript. This work was supported by the Grant Agency of the Czech Republic (grant no. 206/08/0890); by the Internal Grant Agency from Palacký University (IGA_PrF_2014001, IGA_PrF_2015_001, IGA_PrF_2016_001); and by the National Program of Sustainability I, MEYS (grant no. LO1204 – Sustainable development of research in the Centre of the Region Haná).