Completing the hybridization triangle: the inheritance of genetic incompatibilities during homoploid hybrid speciation in ragworts (Senecio).

A new homoploid hybrid lineage needs to establish a degree of reproductive isolation from its parent species if it is to persist as an independent entity, but the role hybridization plays in this process is known in only a handful of cases. The homoploid hybrid ragwort species, Senecio squalidus (Oxford ragwort), originated following the introduction of hybrid plants to the UK approximately 320 years ago. The source of the hybrid plants was from a naturally occurring hybrid zone between S. aethnensis and S. chrysanthemifolius on Mount Etna, Sicily. Previous studies of the parent species found evidence for multiple incompatibility loci causing transmission ratio distortion of genetic markers in their hybrid progeny. This study closes the hybridization triangle by reporting a genetic mapping analysis of the remaining two paired cross combinations between S. squalidus and its parents. Genetic maps produced from F2 mapping families were generally collinear but with half of the linkage groups showing evidence of genomic reorganization between genetic maps. The new maps produced from crosses between S. squalidus and each parent showed multiple incompatibility loci distributed across the genome, some of which co-locate with previously reported incompatibility loci between the parents. These findings suggest that this young homoploid hybrid species has inherited a unique combination of genomic rearrangements and incompatibilities from its parents that contribute to its reproductive isolation.

Introduction

Hybridization is an important contributor to biodiversity and speciation with approximately 25 % of all plant species and 10 % of all animal species estimated to have experienced hybridization during their evolution (Mallet 2005; Baack and Reiseberg 2007). Beyond introgression of novel genetic diversity, the contribution of hybridization to speciation is of particular interest. Genome duplication following hybridization (allopolyploidy) is a frequent mode of speciation in plants with increasing evidence for its occurrence in animals (Rieseberg and Willis 2007; Mable ; Soltis ). Despite the challenges involved in identifying and confirming cases of hybrid speciation without change in ploidy, i.e. homoploid hybrid speciation, improving genetic technologies are accelerating the rate of identification of examples of homoploid hybrid species (Abbott ; Yakimowski and Rieseberg 2014; Goulet ). Moreover, textbook examples of homoploid hybrid speciation are now available for Heliconius butterflies (Mavarez ; Duenez-Guzman ) and Helianthus sunflowers (Rieseberg 2001; Burke ; Lai ), although proof of homoploid hybrid speciation remains controversial (Schumer ; Nieto Feliner ).

Homoploid hybrid speciation is theoretically challenging to explain, particularly when hybrids show sympatry with their progenitors, because reproductive barriers will usually be incomplete, and ongoing hybridization and gene flow are expected to limit the evolution of reproductive isolation and the origin of a homoploid hybrid species (Abbott ; Schumer ). However, various evolutionary processes have been identified as contributing to hybrid speciation without ploidal change. Hybridization is effective at generating a range of new trait combinations and transgressive trait expression that occasionally enable hybrids to exhibit higher fitness than parents in particular ecological contexts (Buerkle ; Lexer ; Schwarz ; Jiggins ; Stelkens and Seehausen 2009). Under these conditions, positive selection can promote the establishment and persistence of new homoploid hybrid species even in the presence of ongoing gene flow (Buerkle ). This process is further facilitated if novel hybrid traits, such as changes in reproductive phenology, pollinator or mating preference, directly reduce gene flow between parents and their hybrids (Servedio ; Marques ).

Another mechanism reducing gene flow and promoting reproductive isolation between hybrids and their progenitors concerns genetic incompatibilities caused by chromosomal rearrangements and/or negative epistasis between parental alleles in the hybrid. New hybrids will often show reduced fitness due to (1) meiotic problems caused by possession of different parental chromosomal rearrangements and (2) the negative interaction of new combinations of alleles inherited from both parents. However, fitness can be recovered through re-assortment of both chromosomal rearrangements and genetic incompatibilities in later generation hybrids (Grant 1981; McCarthy ; Lai ; Schumer ).

Although genetic incompatibilities might be of individually small effect, they appear to be common and their combined effects can be potent evolutionary drivers of speciation (Orr and Presgraves 2000; Presgraves 2010; Schumer ; Lindtke and Buerkle 2015). For example, alleles at different loci that have negative epistatic interactions (Bateson–Dobzhansky–Muller incompatibilities or BDMs) arise readily between isolated populations either through divergent selection or through genetic drift (Coyne and Orr 2004; Paixão ; Schumer ). Recent theoretical models further suggest that selection against negative epistasis within hybrid populations can lead to fixation of combinations of alleles that show genetic incompatibility between hybrids and both parents, thus promoting reproductive isolation (Paixão ; Schumer ). In reality, multiple interacting evolutionary processes probably interact to promote homoploid hybrid speciation. For example, hybrid Helianthus species clearly show contributions from both hybridization-induced genomic reorganization and adaptation to novel ecological niches (Rieseberg 2001; Lexer ; Burke ; Lai ).

Genetic mapping studies are effective at giving a genome-wide perspective of the evolutionary processes driving homoploid hybrid speciation. This study returns to the recent homoploid hybrid origin of S. squalidus (Oxford ragwort) from its parental species, S. aethnensis and S. chrysanthemifolius. The human-aided translocation of hybrids to the UK from a natural hybrid zone on Mount Etna, Sicily, approximately 320 years ago was crucial to the origin of S. squalidus as a new stabilized hybrid species, and its subsequent invasive spread (Abbott ). While geographical isolation allowed the establishment of this homoploid hybrid species, the contribution of genetic incompatibilities to reproductive isolation at the early stages of hybrid speciation is still of interest.

Senecio squalidus shows molecular genetic and quantitative trait divergence from each parent species (and hybrids occurring on Mount Etna), as well as local adaptation associated with latitude within the UK and between the UK and Sicily (Allan and Pannell 2009; Brennan ; Ross 2010). The hybrid zone on Mount Etna is stable despite relatively high gene flow between parent species since they first diverged during the last 150k years (Filatov ) due to both intrinsic hybrid incompatibilities and strong divergent ecological selection associated with elevation (Brennan ; Ross ; Chapman , 2016; Brennan ). All three species are readily hybridized in the glasshouse with few apparent fertility or fitness consequences apart from low seed germination at the F3 generation (Hegarty ; Brennan ). However, recent genetic mapping studies using F2 families derived from crosses between S. aethnensis and S. chrysanthemifolius have characterized genetic incompatibilities in the form of transmission ratio distortion (TRD), breakdown of fitness at the F2 generation, and associations between transmission ratio distortion loci (TRDL) and quantitative traits (Chapman , 2016; Brennan , 2016). These characteristics of hybrid crosses function as genetic incompatibilities by increasing reproductive barriers between taxa (Hall and Willis 2005; Moyle and Graham 2006). Here, we conduct, for the first time, genetic mapping studies using F2 crosses between each of these two species and their hybrid descendent, S. squalidus. In contrast to the hybrid zone on Mount Etna, where selection against hybrids prevents individual hybrid lineages from persisting for many generations (Brennan ), S. squalidus in the UK is a stabilized hybrid descendent that has adapted to a new environment (Abbott ). This set of hybridizing taxa is therefore of considerable interest for a better understanding of the early stages of hybrid speciation.

The genetic mapping studies reported here complete the hybridization triangle of all paired cross combinations between the two parents and their hybrid derivative, and thus increase our understanding of hybrid evolution in this system. Specifically, we test the following hypotheses based on prior knowledge about this hybridizing Senecio system and predictions from the hybrid speciation literature. (1) Hybrid speciation has been accompanied by little major genomic restructuring as suggested by similar genome structures of the parent species; (2) intrinsic genetic incompatibilities are present between S. squalidus and its parents; and (3) they are likely to have been inherited in S. squalidus from its parent species. We discuss our findings in terms of their wider implications for understanding hybridization and homoploid hybrid speciation.

Material and Methods

Mapping families

F2 mapping families were founded from each of three paired cross combinations between three glasshouse grown (F0) individuals representing each of S. aethnensis, S. chrysanthemifolius, and S. squalidus. Senecio aethnensis and S. chrysanthemifolius F0 parental individuals were originally sampled as seed from populations VB and C1 on Mount Etna as described in James and Abbott (2005), whereas the S. squalidus F0 individual was sampled as seed from the Oxford (Ox), UK, population as described in Hiscock (2000). Reciprocal controlled crosses were performed between parental individuals by gently brushing together open flower heads and excluding illegitimate pollen transfer with pollination bags before and after pollination as described in Hiscock (2000) to produce F1 families where the maternal and paternal species of each individual were known. Floral emasculations have been shown to be unnecessary for these typically strongly self-incompatible species (Hiscock 2000). Seeds resulting from these crosses were grown to flowering stage and further reciprocal crosses were performed between full-sib F1 individuals with distinct maternal cytoplasms (i.e. F1 progeny of the same parental individuals but produced from opposite cross directions). From each of the three originally paired species combinations, one reciprocally crossed pair of F1’s was chosen to found each of three F2 mapping families, maintaining approximately equal frequencies of maternal cytoplasm per family. The family hereafter referred to as F2AC was derived from the original cross between S. aethnensis and S. chrysanthemifolius, F2AS was derived from the cross between S. aethnensis and S. squalidus, and F2CS was derived from the cross between S. chrysanthemifolius and S. squalidus. The F2 individuals were labelled in a way to keep track of their maternal origin. The analysis of the F2AC family was previously described in Brennan but is included here for completeness and comparison with the analysis of the two other F2 families.

Genotyping

Genomic DNA was extracted from leaf tissue samples as described in Brennan from all F0, F1 and F2 plants. The number of F2 offspring from which DNA was extracted was 100, 100 and 107 for the F2AC, F2AS and F2CS mapping families, respectively. Approximately 10 % of samples were extracted twice to serve as quality controls.

Samples were genotyped for eight selective primer combinations of Amplified Fragment Length Polymorphisms (AFLPs) according to Brennan , with the final choice of primer combinations and bands scored based on polymorphism in the F0 and F1 parents, high scorability (fluorescence signal >100 rfu), and high repeatability (repeated samples >95 % similar). In addition, the F2 mapping families were genotyped for a total of 75 codominant genetic markers that were found to be polymorphic in the F0 and F1 parents. These comprised 61 expressed sequence tag simple sequence repeats (EST SSRs) and EST indels that were developed from the Senecio expressed sequence tag database (www.seneciodb.org, Hegarty ), and 14 other codominant SSRs and indels derived from previously published Senecio sequences as described in Brennan . Genetic markers were amplified using a three primer system with universal M13 primers fluorescently labelled with FAM6, HEX or NED using a common Polymerase Chain Reaction (PCR) protocol for genotyping on a Beckman Coulter CEQ 8000 capillary sequencer as described in Brennan .

Genetic mapping

Genetic maps were constructed for each of the mapping families using JoinMap v4.0 (Van Ooijen 2001) as previously described for the F2AC mapping family (Brennan ). Because each F2 mapping family was derived from two outbred F0 parents, between two and four alleles segregated at each polymorphic genetic locus. Genotype data for AFLPs and other genetic markers were first formatted according to JoinMap outcrossed mapping family (Crossed Parent or CP type) that allows markers with different segregation patterns, phases and dominance expression to be combined for genetic map construction. Linkage groups (LGs) were identified as sets of markers sharing at least one logarithm of odds (LOD) linkage score >3 and genetic distances <20 Kosambi centiMorgans (cM) following LG regression mapping. Marker order for each LG was determined by iterative rounds of regression mapping, excluding markers that had a large influence on marker order or goodness-of-fit statistics. Each genetic map was summarized for total length according to two different estimators (Chakravarti , Fishman ), map coverage according to Fishman , and genetic marker clustering using dispersion tests (Brennan ).

Genomic rearrangements

Equivalent LGs in each of the three genetic maps were identified on the basis of shared genetic markers. A few genetic markers were found to be mapped to non-equivalent LGs, possibly representing genomic translocations. These markers were excluded from subsequent analyses of relative marker order (synteny). The relative orientation of equivalent LGs on different genetic maps with five or more shared markers was identified by comparing the results of Kendall’s tau correlation tests of shared genetic marker rank order in R v3.1.2 (R Development Core Team 2017). The relative orientation of LGs with fewer than five shared markers was determined by comparing the mean absolute difference in marker rank order. Overall synteny of genetic maps was then assessed using Kendall’s tau correlation tests. Potential genomic rearrangements were identified by examining marker order differences for each paired LG comparison where no combination of LG orientations across the three maps could counteract uncorrelated marker orders in at least one map. Rearrangments typically manifested as transversions with the sequence of marker order differences switching from positive to negative or vice versa. The start and end points of rearrangements were estimated as half the distance between the outermost rearranged marker and the next non-rearranged marker. Recombination rates were compared between rearranged and non-rearranged genomic regions in each map using two sample unpaired rank-sum Wilcox tests of mapped marker distances.

Transmission ratio distortion

To identify genomic regions containing genetic incompatibilities between species, genotype data for each marker in each mapping family were initially tested (Test 1) for TRD against the null hypothesis of Mendelian segregation with χ2 tests using Microsoft Excel 2010 (Microsoft Corp). Genetic markers were considered to show TRD if χ2 tests were significant at a 95 % confidence level before examining the effect of using a 100-fold stricter confidence level to account for multiple marker tests per map.

All markers linked to a genetic incompatibility locus are expected to show similar patterns of TRD. Therefore, clusters of genetic markers showing similar patterns of TRD and linked by less than 10 cM map distance were identified as TRDLs. Single genetic markers showing independent patterns of TRD and located more than 10 cM from other markers with TRD were considered as additional TRDLs, albeit with less supporting evidence. The extent of a genomic region affected by TRD was estimated as halfway to the next marker not showing TRD.

Genotype data were further examined to determine the possible causes of TRD by means of the following tests. Test 2: Cytonuclear incompatibilities dependent on cross direction were examined by testing for genotypic TRD in the subsets of each mapping family in each reciprocal cross direction (i.e. sharing the same species cytoplasm). Test 3: For those TRDLs containing genetic markers for which all alleles could be assigned to each F0 parent, the influence of pre-zygotic incompatibilities acting at the haploid gametophyte stage was assessed by performing χ2 tests of allelic frequencies against null expectations. Dominantly scored markers or markers where parental genotypes shared some alleles could not be tested in this way. Test 4: Also, for those TRDLs containing genetic markers for which all alleles could be assigned to each F0 parent, deficits or excesses of heterozygotes in terms of F0 parental origin were tested with χ2 tests against null expectations. Test 5: Two-locus negative epistatic interactions, also known as BDMs were tested by building contingency tables of genotype combinations for the most distorted marker within each pair of TRDLs and performing Fisher’s exact tests for biases in genotype frequency combinations.

Co-location of TRDLs across the different genetic maps and with genomic rearrangements

The locations of observed TRDLs and genomic rearrangements were compared between the different genetic maps by transposing them onto the F2AC map according to common genetic markers. For this, the location of the genetic marker in the F2AC map that showed most TRD in each TRDL and the midpoint of each rearranged genomic region were used for subsequent analyses of co-location. Co-location between TRDLs identified by genotype tests in each pair of genetic maps after transposition onto the F2AC genetic map was tested using sampling without replacement tests as described in Brennan . Briefly, the genetic map was divided into ‘n’ intervals of equal length and the frequency of TRDL occurrence and co-occurrence in each interval was tested against null expectations of no association between the genomic distributions of TRDLs in each genetic map. The influence of including either TRDs with individual- and map-level 95 % confidence levels at a range of interval sizes on the test statistic was investigated. Those TRDLs that mapped to different equivalent LGs when transposed onto the F2AC genetic map were excluded from the TRDL co-location analysis due to uncertainty over their map position. Following these analyses, the co-location between genomic rearrangements and TRDLs was examined in the same way.

Results

Genetic maps

Details of all genetic markers used to genotype individuals for the construction of genetic maps including: type (SSR, indel or AFLP), EST match in the Senecio EST database, primer sequences, F0 genotypes and genetic mapping information are provided in . The full F2AS, F2AC and F2CS genetic maps are illustrated in with a summary of five LGs presented in Fig. 1 and summary statistics provided in Table 1. The F2AC genetic map was described previously (see Brennan and Dryad Digital Repository: doi:10.5061/dryad.7b56k) and is included here for comparative purposes. F2 genotype data used to construct the F2AS and F2CS genetic maps are available at (DRYAD Digital Repository; doi:10.5061/dryad.82d5f33).

Figure 1.

Interleaved genetic maps of selected linkage groups (LGs) from the F2AC, F2AS and F2CS mapping families. Map distances in Kosambi centiMorgans are shown in the scale to the left of LGs. Linkage groups are represented by vertical bars with mapped marker positions indicated with horizontal lines. Linkage group names are presented in bold above each LG with letters indicating the pair of F0 species; S. aethnensis (A), S. chrysanthemifolius (C) or S. squalidus (S) that founded the mapping family and numbers indicating equivalent LGs that share genetic markers across the three maps. Weakly linked LGs that are thought to belong to the same chromosome are aligned vertically under a single overall LG name. Marker names are listed to the left of LGs in grey if they are common to another genetic map or in black if they are uniquely present on that genetic map. Dotted lines link common marker positions on the equivalent LGs of different genetic maps. Black shaded portions of LGs indicate chromosomal transversions identified from switches in marker order compared with equivalent LGs. See for a depiction of all LGs corresponding to 10 chromosomes.

Table 1.

Summary genetic linkage map statistics for the F2AC, F2AS and F2CS genetic maps

Linkage group	Length			No. of genetic markers			No. of codominant markers			No. dominant markers			Add 2S length			Method 4 length
	F2AC	F2AS	F2CS	F2AC	F2AS	F2CS	F2AC	F2AS	F2CS	F2AC	F2AS	F2CS	F2AC	F2AS	F2CS	F2AC	F2AS	F2CS
1	44.5	47.0	46.4	18	18	18	14	11	12	4	7	6	50.1	51.8	51.0	49.7	52.5	51.9
2	29.1	46.5	39.1	7	14	12	3	4	5	4	10	7	34.7	51.3	43.7	38.8	53.7	46.2
3	42.6	21.7	25.3	10	10	15	4	3	5	6	7	10	48.2	26.5	29.9	52.1	26.5	28.9
4A	41.3	39.9	67.2	10	14	12	4	5	7	6	9	5	46.9	44.7	71.8	50.5	46.0	79.4
4B	–	6.6	–	–	4	–	–	2	–	–	2	–	–	11.4	–	–	11.0	–
5A	25.8	26.9	23.0	9	7	12	4	3	5	5	4	7	15.1	31.7	27.6	15.8	35.9	27.2
5B	9.5	13.2	6.7	4	3	2	4	3	2	0	0	0	47.3	18	11.3	52.1	26.4	20.1
6A	41.7	28.3	21.0	9	16	5	2	5	0	7	11	5	19.8	33.1	25.6	16.6	32.1	31.5
6B	–	–	8.9	–	–	8	–	–	5	–	–	3	–	–	13.5	–	–	11.4
7A	14.2	12.3	15.9	13	10	15	1	2	1	12	8	14	19.8	17.1	20.5	16.6	15.0	18.2
7B	3.2	–	3.2	2	–	2	1	–	1	1	–	1	8.8	–	7.8	9.6	–	9.6
8A	27.5	24.1	15.7	22	18	12	6	7	2	16	11	10	33.1	28.9	20.3	30.1	26.9	18.6
8B	5.2	–	5.9	2	–	2	1	–	0	1	–	2	10.8	–	10.5	15.6	–	17.7
9	15.0	7.1	5.4	8	8	12	1	1	1	7	7	11	20.6	11.9	10.0	19.3	9.1	6.4
10A	10.0	11.4	9.2	11	11	14	3	3	4	8	8	10	15.6	16.2	13.8	12.0	13.7	10.6
10B	4.2	4.3	1.9	2	4	2	2	3	2	0	1	0	9.8	9.1	6.5	12.6	7.2	5.7
total	313.8	289.3	294.8	127	137	143	50	52	52	77	85	91	391.6	351.2	364.3	407.1	356.0	383.3
mean	22.41	22.25	19.65	9.07	10.54	9.53	3.57	4.00	3.47	5.50	6.54	6.07	27.97	27.01	24.28	29.08	27.39	25.56
stdev	15.56	14.91	18.52	5.88	5.22	5.55	3.37	2.61	3.20	4.55	3.69	4.33	15.56	14.91	18.52	16.65	16.12	20.18
unmapped unlinked				9	6	4	1	2	2	8	4	2
unmapped problematic				9	8	9	3	2	1	6	6	8
unmapped total				18	14	13	4	4	3	14	10	10

Linkage groups are named as in Brennan where ‘A’ and ‘B’ after numbers indicate that these LGs probably belong to the same chromosome. ‘–’ indicates that an equivalent linkage was not observed in a genetic map; this is sometimes due to mapping of these markers to the preceding LG. Map length measures are in Kosambi centiMorgan units. Add 2S length is an estimate of chromosome length calculated as LG length plus twice mean LG distance. Method 4 length is another estimate of chromosome length calculated as LG length times (marker number + 1)/(marker number − 1).

Each map comprised 13 to 15 LGs based on our minimum linkage criteria of greater than 3 LOD linkage between pairs of markers and <20 Kosambi cM distance between markers. However, some pairs of LGs showed weaker evidence of linkage or were part of the same equivalent LG in another genetic map leading to the conclusion that these maps correspond to the n = 10 chromosomes expected for these Senecio species (Alexander 1979). The F2AS and F2CS maps contained slightly more markers (139 and 143, respectively) than the F2AC map (127 markers), but were slightly shorter in total map length (289.3 and 294.8 Kosambi cM, respectively) relative to the F2AC map (313.8 Kosambi cM). Genetic markers were separated by mean distances of 2.8, 2.4 and 2.3 Kosambi cM within the F2AC, F2AS and F2CS genetic maps, respectively, with between >95.6 % and > 99.8 % of the genome predicted to be within 5 and 10 Kosambi cM of a mapped maker. However, according to dispersion tests the distributions of genetic markers were significantly clumped across all three maps (P < 10–16, all maps), indicating that some regions of each map show better marker coverage than others.

Paired comparisons between the three genetic maps showed that almost half (62 to 70) of the genetic markers were shared between each pair of maps allowing identification and orientation of equivalent LGs [see and Table 2]. The maps showed high synteny overall as indicated by high Kendall’s tau correlation coefficients for genetic marker rank order (Fig. 2). Based on the magnitude of the tau coefficient, the overall AS-CS comparison showed highest synteny, followed by AC-AS and then AC-CS (Fig. 2). Nonetheless, one to three genetic markers per map were discordant with the overall sharing of markers between LGs and were present on a non-equivalent LG in another of the maps, as labelled in Fig. 2. Moreover, ten instances of switches in marker order, corresponding to five genomic regions, were detected when shared marker order was examined between pairs of maps at an LG level (Fig. 1, Table 2). These likely genomic rearrangements were present on LGs 8 and 4 in the F2AC and F2AS maps, respectively, and on LGs 3, 9 and 10 in the F2CS map, as highlighted in Fig. 1. It was further evident that mapped marker distances were significantly shorter within genomic regions associated with these rearrangements compared with other genomic regions (mean ± SD 0.6 ± 0.7 versus 3.8 ± 4.1 for F2AC, 0.6 ± 0.8 versus 3.1 ± 3.8 for F2AS, 0.5 ± 0.6 versus 3.0 ± 4.3 for F2CS, Wilcox test P < 10–05 for each map) in accordance with expectations of reduced recombination within rearranged regions.

Table 2.

Summary of transversions detected by tests for synteny of shared marker order for each LG in each paired genetic map comparison.

Map comparison	# markers	Tau (P value)	absdiff- abs(diff)	AC LG	AC start marker	AC start position (cM)	AC end marker	AC end position (cM)	length (cM)	transverted map	min start position (cM)	max end position (cM)	length (cM)
F2AC and F2CS	6	0.69 (0.06)	0.67	AC3	E5M3_405	36.5	E1M5_202	41.3	4.8	F2CS	34.2	41.3	7.1
F2AS and F2CS	8	0.21 (0.55)	2.25	AC3	E5M6_187	34.2	E5M6_453	36.4	2.2
F2AS and F2CS	7	0.1 (0.76)	2.29	AC4	E4M7_152	17.3	E4M7_179	20.6	3.3	F2AS	16.1	20.6	4.5
F2AS and F2CS	8	0.57 (0.06)	1	AC4	E8M5_157	16.1	EC512	18.1	2
F2AS and F2CS	11	0.38 (0.12)	2.36	AC8	E5M6_401	13.6	EC482	17.5	3.9	F2AC	13.6	17.5	3.9
F2AC and F2CS	9	0.56 (0.04)	1.33	AC8	E1M3_121	14.8	E5M6_103	16.6	1.8
F2AC and F2CS	6	0.33 (0.47)	1.33	AC9	V45	5.8	E1M3_317	12.7	6.9	F2CS	5.8	15.7	9.9
F2AS and F2CS	6	0.46 (0.22)	1	AC9	E8M5_153	12.3	E8M5_205	15.7	3.4
F2AC and F2CS	8	0.49 (0.10)	1.25	AC10	E1M7_211	2.8	ES20	4.4	1.6	F2CS	2.8	5.4	2.6
F2AS and F2CS	9	0.43 (0.11)	1.56	AC10	E8M7_223	5.1	EC688	5.4	0.3

# markers is the number of common markers on the equivalent LG of each pair of compared genetic maps. tau is the Kendall paired rank correlation test summary statistic, values range from 0 to 1 with larger values indicating higher synteny. P value is the probability of the observed tau values, test results greater than 0.05 are shown indicating that the compared marker orders are insignificantly different than random. absdiff-abs(diff) summarizes marker order differences, it is the mean marker absolute difference in rank order minus the absolute of mean marker differences in rank order, larger absdiff-abs(diff) values indicate transversions in marker order that start with negative differences and end with positive difference that cancel each other out leading to smaller abs(diff) values. AC LG is the equivalent F2AC LG where the difference in marker order was found. AC start/end marker is the nearest F2AC mapped marker to the start/end of the change in marker order. AC start/end position is the equivalent F2AC map start/end position of the change in marker order, when start/end markers are not present on the F2AC map; an approximate position is calculated as the distance between the start/end marker to the nearest F2AC mapped marker. Length (cM) is the centiMorgan distance between the F2AC map start and end of the change in marker order. Transverted map is the genetic map that shows marker order differences in the same map region with both other maps. Min start position and max end position are the smaller of two tranversion start positions and the larger of the two tranversion end positions, respectively, in the same genetic map region based on the two paired map comparisons.

Figure 2.

Paired comparisons of marker order between the F2AC, F2AS and F2CS genetic maps. Points indicate the relative map order of genetic markers common to each pair of compared genetic maps. Blocks of alternating light grey and black points indicate different LGs. Labelled circles indicate genetic markers in discrepant map positions (≥10 difference in relative map order) in each pair of compared genetic maps. The dashed lines indicate identical marker order for comparison. Summaries of Kendall’s rank correlation tests are shown in the top left corner of each panel.

Results of all TRD tests for genetic markers analysed in the F2AS and F2CS families are presented in and , respectively. Equivalent results for the F2AC family are available in Table S2 of Brennan . The three F2 mapping families differed significantly in the frequency of genetic markers showing genotypic TRD [Test 1, ], with 26.8 %, 14 % and 2.9 % of markers exhibiting TRD in the F2AC, F2CS and F2AS mapping families, respectively. When genotypic TRD was tested using the more stringent map-wide 95 % confidence level, no difference was detected between the F2AC and F2CS mapping families in frequency of genetic markers showing TRD, but the difference between these two families and the F2AS mapping family remained, with the latter map containing no markers exhibiting genotypic TRD at this threshold [see]. Unmapped genetic markers were more likely to show genotypic TRD in each mapping population.

After combined markers with genotypic TRD into TRDLs, it was apparent that fewer TRDLs were detected in the F2AS and F2CS maps (containing four and six TRDLs, respectively) (Table 3, Fig. 3) than in the F2AC genetic map (containing nine TRDLs) (Table 3 in Brennan ). Although most TRDLs identified in the F2AS and F2CS maps mapped to equivalent LGs in the F2AC genetic map, there were two notable exceptions. These included a TRDL represented by genetic marker E3M5_65, which mapped to LGs AC2 and AS4, and another represented by E8M5_110 that mapped to LGs AC8 and CS1. These exceptions could represent translocations or errors in mapping potentially caused by TRD itself.

Table 3.

Summary of TRDLs observed for tests of different TRD mechanisms in the F2AS and F2CS genetic maps

Linkage group_map position (cM)	Equivalent AC LG and map position (cM)	Genetic marker with greatest TRD	Genotypic TRD	Asymmetric / cytonuclear TRD	Pre-zygotic TRD	Heterozygote deficit	Epistasis / BDMs (minority genotype)
AS2_3.7 singleton	AC2_9.3	EC978	Yes	No	No	Yes	No
AS4A_37.9 singleton	AC2_11.3	E5M3_65	Yes	Squal	–	–	No
AS3_0.0 cluster	AC3_29.5	ES1	Yes	Squal	No	Yes	No
AS5A_26.9 singleton	AC5A_12.6	EC1470	yes	Squal	No	–	No
CS1_6.6 cluster	AC1_0.0	EC296B	Yes	No	Squal	Yes	CS1_6.7 BC, CS4A_38.1 BD
CS4B_0.0 singleton	AC4A proximal	E5M3_219	Yes	No	–	–	No
CS4A_38.1 singleton	AC4A_41.3	E1M3_264	Yes	No	–	–	CS1_6.6 DB, CS1_6.7 CD
CS6A_0.0 singleton	AC6A proximal	E5M6_397	Yes	No	Yes	–	No
CS7A_unmapped	AC7A_1.0	E1M5_269	Yes	No	Yes	–	No
CS1_6.7 cluster	AC8A_27.5	E8M5_110	Yes	No	Yes	–	CS1_6.6 CB, CS4A_38.1 DC
AS1_28.4 singleton	AC1 central	E5M3_104	No	Aeth	–	–	–
AS2_29.7 singleton	AC2 central	ES56	No	Aeth	No	Yes	–
AS3_16.9 singleton	AC3 distal	ES58	No	Aeth	No	Yes	–
AS6_6.0 singleton	AC6 proximal	E1M5_140	No	Aeth	No	–	–
AS1_37.8 cluster	AC1_27.1	EC74	No	Squal	Squal	No	–
AS2_44.6 singleton	AC2 distal	ES74B	No	Squal	No	–	–
AS5A_0.0 cluster	AC5A_3.5	E1M3_254	No	Squal	–	–	–
AS10A_7.6 cluster	AC10A central	E8M7_223	No	Squal	No	–	–
CS4B_11.5 cluster	AC4A_18.9	ES2	No	Chrys	No	Yes	–
CS4B_67.2 singleton	AC4A distal	E1M8_106	No	Chrys	–	–	–
CS7A_13.5 singleton	AC7A distal	E8M7_226	No	Chrys	No	–	–
CS8A_6.8 singleton	AC8A central	E1M5_88	No	Both	–	–	–
CS10A_9.2 cluster	AC10A_4.4	ES20	No	Chrys	no	Hom	–
CS2_2.0 singleton	AC2 proximal	E1M3_163	No	Squal	–	–	–
CS3_21.7 cluster	AC3 distal	E5M6_187	No	Squal	No	–	–

Linkage group names correspond to Fig. 1. All tests for TRD involved χ2 tests against null expectations at a per-marker 95 % confidence level unless stated otherwise. ‘-’ indicates that genetic marker genotypes did not allow particular TRD tests. Cluster indicates that the TRDL is represented by multiple linked genetic markers less than 10 cM apart, whereas singleton indicates that the TRDL is represented by a single marker. Equivalent F2AC genetic map positions are estimated as ‘proximal’, ‘central’ or ‘distal’ based on additional linked genetic markers when the marker with greatest TRD was not itself present on the F2AC map. Genotypic TRD indicates markers that showed TRD for genotype frequencies. Asymmetric/cytonuclear TRD indicates if TRD was present for one parental cytotype only, showing the affected cytotype as ‘Aeth’ for S. aethnensis, ‘Chrys’ for S. chrysanthemifolius and ‘Squal’ for S. squalidus. Pre-zygotic TRD indicates if TRD was present for allelic frequencies showing the minority parental allele for significant cases. Heterozygote deficit indicates if a significant deficit of heterozygotes with both parental alleles was observed, with ‘Hom’ indicating the opposite case of a significant deficit of homozygotes. Epistasis / BDMs indicates if Fisher’s exact tests of paired genetic marker genotype combinations showed significant interactions indicative of two-locus BDM incompatibilities. Epistatic interactions are summarized as the other interacting TRDL locations and the minority genotype at the ‘home’ and ‘away’ TRDL. Genotypes are described in mapmaker format where ‘A’ and ‘B’ indicate homozygous parental alleles (parents listed alphabetically), ‘H’ indicates heterozygous parental alleles, ‘C’ indicates not homozygous A parent and ‘D’ indicates not homozygous B parent.

Figure 3.

Genetic map locations and evidence for co-location of F2AS and F2CS transmission distortion loci transposed onto the F2AC genetic map. The F2AC genetic map was used as a reference against which to compare TRDL locations. Genetic map features are as described in the notes to Fig. 1. TRDL locations are represented as vertical lines to the right of F2AC LGs and named after the genetic map in which they were observed. The lines depicting TRDLs extend to cover linked genetic markers showing TRD. A bold cross hatch indicates the location of the marker with greatest TRD within that TRDL. Black lines indicate TRDLs identified from significantly biased genotype frequencies while grey lines indicate TRDLs identified by other tests (see Table 2). * after TRDL names indicates that χ2 test resulted in P < 0.0005. ~ before TRDL names indicates that the marker with greatest TRD was not found on the F2AC map but the approximate position and extent of the TRDL was estimated from synteny across genetic maps. ^ before TRDL names indicates that the marker with greatest TRD was present on a different LG compared to its F2AC location (see Table 2).

Further tests indicated the possible causes of TRD at a locus. Thus, Test 2 indicated examples of asymmetric TRD (i.e. dependent on cross direction) at three of four TRDLs resolved in the F2AS mapping family, but at none of the six loci resolved in the F2CS family (Table 3). Rather surprisingly, some additional asymmetric TRDLs (eight and four in the F2AS and F2CS maps, respectively) were detected by these cross-specific tests (Table 3, Fig. 3) that were not detected in the analysis of entire F2 progenies (Test 1).

Tests of pre-zygotic (or allelic) TRD (Test 3) were significant for one F2AS and four F2CS TRDLs (Table 3), whereas tests of parental allele heterozygote surplus or deficit (Test 4) showed deficits for one F2CS and two F2AS TRDLs resolved in the analyses of entire F2 families and two additional F2AS TRDLs resolved in the tests of asymmetric TRD (Table 3). Finally, tests of negative epistasis between TRDLs (Test 5) were significant for three interacting F2CS TRDLs represented by markers EC296B, E1M3_264 and E8M5_110 (Table 3). The minority genotype combinations for these TRDLs suggested a lack of S. squalidus-like genotypes. Two of these TRDLs, EC296B and E8M5_110, while on separate LGs in the F2CS map, form part of a large TRDL located on the AC1 LG of the F2AC map (Brennan ).

Co-location of TRDLs (genetic incompatibilities) across genetic maps and with genomic rearrangements

Comparison of the genomic distribution of TRDLs indicated that all LGs contained TRDLs across the three genetic maps with between one and five TRDLs observed per LG (Fig. 3, Table 3). There was evidence for co-location of TRDLs between genetic maps. In particular, the F2AC and F2CS genetic maps showed significantly more co-located TRDLs than expected by chance for all subsets of TRDLs and across a range of map interval test lengths from 2 to 8 cM [see]. Co-location tests at larger map interval lengths were probably more conservative as the random chance of multiple TRDLs occurring within the same large interval increased. Visual inspection of the maps identified these co-located TRDLs on LGs 1, 4, 7 and 10. No significant co-location of TRDLs was evident from comparisons of the F2AC and F2AS or the F2AS and F2CS maps. Co-location tests involving TRDLs that remained significant after correction for multiple testing found significant co-location of TRDLs on LGs 1 and 4 of the F2AC and F2CS maps [see. No evidence was found for TRDLs being co-located with the midpoints of genomic rearrangments as inclusion of these data into the co-location analyses did not change the significance of test statistics [see].

Discussion

The genetic mapping and TRD analyses conducted on F2 families generated from pairwise species crosses indicate that the homoploid hybrid species, Senecio squalidus, has inherited genetic incompatibilities from both of its parental species, S. aethnensis and S. chrysanthemifolius. These findings lend support to a model of how reproductive isolation of a homoploid hybrid species can be initiated by inheritance of pre-existing genetic incompatibilities between the parental species (Schumer ). Although reproductive isolation of S. squalidus from its parents is primarily dependent on ecogeographic isolation, with some isolating effects possibly resulting from genetic drift or selection during its origin and establishment in the UK (James and Abbott 2005; Abbott ; Ross 2010; Brennan ), our results indicate that inherited genetic incompatibilities also contribute to the reproductive isolation of this homoploid hybrid species.

Our previous genetic mapping study using the same S. aethnensis and S. chyrsanthemifolius parental individuals as in this study suggested that the large-scale structure of the genomes of both species was similar with no genetic evidence of fusions, fissions or translocations among chromosomes (Brennan ). This finding was further supported by the results of an independent genetic mapping study conducted by Chapman . The present investigation extended these analyses to a comparison of the genomic structure of S. squalidus with its progenitors and indicated that some large-scale genomic restructuring had occurred during the origin of S. squalidus, although the three genetic maps were very similar overall, in terms of number and length of LGs detected. Comparisons between the three genetic maps indicated that the 13 to 15 LGs present in each map could be assigned to the ten chromosomes expected for these Senecio species (Fig. 1 and , Table 1). The maps also showed similar total and individual LG lengths, which corroborates the similar 2C nuclear DNA content measures of 1.57, 1.63 and 1.41 pg recorded for S. aethnensis, S. chrysanthemifolius and S. squalidus, respectively (Coyle and Abbott, unpubl. results). It seems, therefore, that genome size increase caused by retrotransposon proliferation activated by hybridization as reported in Helianthus (Baack ; Ungerer ) and Aegilops (Senerchia ) is probably not a feature of this Senecio system.

Nature of genomic rearrangements

Genetic mapping showed that approximately half of all component genetic markers were shared between each pair of genetic maps and these shared markers indicated high overall synteny between maps based on correlation tests (Fig. 2). A few genetic translocations between LGs, involving only one to three genetic markers that were found on non-equivalent LGs per paired map comparison, were observed (Fig. 2), and might reflect small-scale genomic translocations affecting individual genetic markers at a scale of 2 to 5 cM (one to two times the mean cM distance between mapped markers). However, an alternative explanation is that some of these genetic markers that map to different LGs in different genetic maps also show strong genotypic TRD in one or both genetic maps (E5M3_65 for the F2AC versus F2AS comparison and E1M7_207 and E8M5_110 for the F2AC versus F2CS comparison). Strong TRD could cause false associations between genotypes of unlinked markers leading to errors in map location.

Tests of marker order at the LG level found evidence for genomic transversions affecting at least five different genomic regions; one present on each of the F2AC (LG8) and F2AS (LG4) maps and three on LGs 3, 9 and 10 of the F2CS map (Fig. 1). These results suggest that the genomes of the three parental individuals representing each of the three study species are distinguished by one to three genomic rearrangements. Determining which rearrangement is associated with each species and which rearrangements might have been inherited by S. squalidus from its parental species would require additional genetic mapping studies of within species crosses. These rearranged genomic regions also show significantly shorter mapped marker distances than other genomic regions. This is because recombination is negatively selected within these genomic regions in the progeny of individuals heterozyogous for the rearrangement as it generates large deleterious insertion-deletion mutations (Fishman ).

The genomic rearrangements identified could indicate genomic regions that potentially harbour many genes of functional significance that are protected from interspecific gene flow and upon which selection can act to promote divergent evolution. Genomic inversions have been found to be associated with multiple examples of divergent adaptive evolution in the presence of gene flow, for example, perennial and annual ecotypes of Mimulus guttatus (Twyford and Friedman 2015), and different Mullerian host mimics in the butterfly, Heliconius numata (Joron ). It would be of interest to investigate further the potential contribution that genomic rearrangements between S. aethnensis and S. chrysanthemifolius make towards maintaining ecological differentiation despite gene flow across the Mount Etna hybrid zone. The interspecific genomic rearrangements observed in these genetic maps could also contribute to genetic incompatibility between S. squalidus and both of its parental species.

Genomic reconstruction during the origin of S. squalidus agrees with reports of the same in other homoploid hybrid species, e.g. in Helianthus (Burke ) and in Iris (Tang ; Taylor ). In these cases, genomic reorganization was interpreted as a ‘genomic shock’ response to hybridization (McClintock 1984; Rieseberg 2001; Chen and Ni 2006), enabling stabilization of the hybrid genome, and associated positively with evolutionary divergence of the parental genomes. The parental species of S. squalidus are estimated to have diverged relatively recently, around 150 000 years ago, once suitable habitats for the high altitude species S. aethnensis became available with the rise of the volcano, Mount Etna, in Sicily (Chapman ; Osborne ). Despite their very recent origin, some genomic rearrangements appear to have already emerged between these two species. In combination with findings for other homoploid hybrid species (Tang ; Taylor ), our results suggest that genomic restructuring is a frequent feature of the successful establishment of new hybrid species in combination with ecological and/or spatial divergence from parental species (Buerkle ; Baack and Rieseberg 2007, Karrenberg ).

The presence of TRD among genotyped progenies is indicative of genetic incompatibilities between the parental lines because particular alleles, genotypes, or combinations of these have been selected against in hybrid offspring (Fishman ; Harushima ; Hall and Willis 2005; Moyle and Graham 2006). Comparisons between genetic maps showed that TRD was present at multiple genomic regions across all chromosomes. Most genotype-level TRDLs and the highest proportion of genetic markers exhibiting TRD were evident in the F2AC map (Fig. 3, Table 3 and ). This observation fits with evidence that hybrid incompatibilities have accumulated and been reinforced by divergent ecological selection between the parental species on Mount Etna (Osborne ; Brennan ; Chapman ; Filatov ), in contrast to incompatibilities involving the hybrid species that no longer interacts with its parents due to geographic isolation.

The crossing design used to produce the F2 mapping families involved full-sib F1 crosses raising the possibility that some of the observed TRD could be due to bi-parental inbreeding. These species are self-incompatible and S. squalidus has been shown to suffer from inbreeding depression when selfed (Brennan , 2013). Inbreeding depression would affect patterns of TRD in the form of selection against reconstituted homozygous F0 parental genotypes in the F2 progeny. However, at TRDLs where tests of TRD of heterozygosity of marker alleles of different F0 parental origin could be applied, the majority of significant results were in favour of a deficit of heterozygotes (one out of one tests for the F2AC map, four out of four tests for the F2AS map, two out of three tests for the F2CS map, Table 3). Therefore, most TRDLs in these mapping families appear to be caused by genetic incompatibilities between species rather than inbreeding depression.

Less TRD was evident in the F2 families produced from crosses between each parental species and S. squalidus than between the two parent species. For the families involving S. squalidus, 14 % of markers were distributed across six TRDLs in the F2CS family compared with 2.9 % of markers showing TRD across four TRDLs in the F2AS family (). The extent of genotypic TRD was also more pronounced for many F2CS TRDLs relative to F2AS TRDLs, indicating that S. squalidus inherited a greater number of S. aethnensis-like incompatibility alleles or local rearrangements that preferentially cause genetic incompatibility with the S. chrysanthemifolius parent. Asymmetric backcross incompatibility and directions of introgression have been reported for a number of hybridizing species and can be caused by cytoplasmic incompatibilities between nuclear and chloroplastic genomes (Buerkle and Rieseberg 2001; Martin ; Scascitelli ; Senerchia ; Abbott 2017). In the case of S. squalidus, samples have been found to share the same chloroplast DNA haplotype with both parental species, so the direction of the original hybrid cross is currently uncertain (Abbott ; Comes and Abbott 2001; Simon Hiscock, unpubl. data). The TRD tests that took cytoplasmic identity into account (Test 2) found only two instances of asymmetric TRD in each of the F2AS and F2CS mapping families, suggesting that cytoplasmic incompabilities are minor contributors to the overall hybrid incompatibility observed in this system and supporting the hypothesis that hybridization in both cross directions could contribute to gene flow and hybrid evolution.

The greater prevalence of genetic incompatibilities between S. squalidus and S. chrysanthemifolius does not appear to have biased parental contributions to the hybrid genome of S. squalidus (James and Abbott 2005; Brennan ; Filatov ). Instead, the effect of these genetic incompatibilities on hybridization dynamics would seem to be restricted to smaller genomic regions. Considering all the forms of TRD identified, each cross showed multiple TRDLs distributed across the genome that function in a mixture of cross directions, so that their combined effect would contribute to genetic incompatibility in both cross directions.

Similar to results previously reported for the F2AC mapping family (Brennan ), additional tests of TRD provided evidence for cytonuclear incompatibilities, allelic pre-zygotic incompatibilities, heterozygote (and homozygote) deficit and two-locus epistatic incompatibilities as causes of TRD in both F2AS and F2CS mapping families (Table 3). Moreover, neighbouring TRDLs (<10 cM apart) were identified to exhibit TRD resulting from different causes, as demonstrated for the neighbouring pairs of F2AS TRDLs on LGs 1 and 2 (Fig. 3). It seems likely, therefore, that additional TRDLs would be detected if TRD in hybrid crosses were to be studied at greater resolution with more markers employed. Insufficient genomic resolution might, in part, explain why some TRDLs were observed in only one family out of the three studied. It needs emphasizing that the construction of these genetic maps involved a single representative of each species and therefore represents a snapshot of all the genetic incompatibilities that are present in this system. Genetic maps built from different parents might reveal a slightly different subset of genetic incompatibilities if the alleles causing these incompatiblities have not been fixed in the different species, as noted in hybridizing Mimulus guttatus and M. nasutus (Sweigart ; Martin and Willis 2010). It is less likely, but not inconceiveable, that new genetic incompatibilities such as BDMs could have emerged between the parental species and S. squalidus since it became allopatrically isolated in the UK approximately 320 generations ago.

The results of this study confirm that, in addition to previously identified intrinsic hybrid incompatibilities between S. aethnensis and S. chrysanthemifolius (Brennan , 2014; Chapman ), genetic incompatibilities are also present between the homoploid hybrid species, S. squalidus, and its parental species. We also found support for the hypothesis that these genetic incompatibilities in S. squalidus were inherited from its progenitors by testing for genetic incompatibilities that were shared between the different F2 mapping families. TRDLs were found to be significantly co-located between the F2AC and F2CS maps based on four co-located TRDLs on LGs 1, 4, 7 and 10 (Fig. 3). Taken overall, the evidence we obtained of multiple shared TRDLs in the genetic maps of crosses between the hybrid species and its progenitors supports the evolutionary potential of inheritance and re-assortment of hybrid incompatibilities (Grant 1981; Paixão ; Schumer ).

Conclusions

While evidence for the contribution of hybridization to speciation continues to accumulate in the literature (Abbott 2017 and references therein), our understanding of the process at a genomic level is still very limited. Our study addresses this knowledge gap using a genetic mapping approach to investigate the structure of the genome of a new homoploid hybrid species in comparison to its progenitors. Our results reinforce the view that hybridization has heterogeneous effects across the genome at multiple dispersed genomic locations. A challenge for the future is to examine a greater variety of naturally hybridizing systems at a sufficiently dense genomic resolution to determine the generality of these observations and to zoom in on the particular genes or genomic structures acting as hybridizing barriers (e.g. Christe ). There continues to be a need to integrate new genetic data with data on the effects of hybridization on quantitative traits and fitness, particularly in the environments where hybridization actually occurs (Goulet ). The developing applicability of high-throughput sequencing methods and their analysis to non-model hybridizing systems will contribute to these issues and provide new insights into the evolutionary consequences of hybridization.

Data archiving

Mapping family genotype data and genetic map information has been deposited with the DRYAD Digital Repository (https://doi.org/10.5061/dryad.82d5f33). Other results can be found in the supporting information.

Source of Funding

The research was funded by an Natural Environment Research Council (grant NE/D014166/1) to RJA as Principal Investigator.

Contributions by the Authors

A.C.B., S.J.H. and R.J.A. conceived the study. A.C.B. performed the experiments and analysed the results. A.C.B., S.J.H. and R.J.A. contributed to writing the paper.

Conflict of Interest

None declared.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.