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

Abstract 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 et al. 2011;Soltis et al. 2014). Despite the challenges involved in identifying and confirming *Corresponding author's e-mail address: a.c.brennan@durham.ac.uk 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 et al. 2013;Yakimowski and Rieseberg 2014;Goulet et al. 2017). Moreover, textbook examples of homoploid hybrid speciation are now available for Heliconius butterflies (Mavarez et al. 2006;Duenez-Guzman et al. 2009) and Helianthus sunflowers (Rieseberg 2001;Burke et al. 2004;Lai et al. 2005), although proof of homoploid hybrid speciation remains controversial (Schumer et al. 2014b;Nieto Feliner et al. 2017).
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 et al. 2013;Schumer et al. 2014b). 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 et al. 2000;Lexer et al. 2003;Schwarz et al. 2005;Jiggins et al. 2008;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 et al. 2000). 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 et al. 2011;Marques et al. 2016).
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 et al. 1995;Lai et al. 2005;Schumer et al. 2015).
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 et al. 2014a;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 et al. 2014;Schumer et al. 2015). 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 et al. 2014;Schumer et al. 2015). 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 et al. 2003;Burke et al. 2004;Lai et al. 2005).
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 ). 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 et al. 2012;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  due to both intrinsic hybrid incompatibilities and strong divergent ecological selection associated with elevation Ross et al. 2012;Chapman et al. 2013Chapman et al. , 2016Brennan et al. 2014). All three species are readily hybridized in the glasshouse with few apparent fertility or fitness consequences apart from low seed germination at the F 3 generation (Hegarty et al. 2008;Brennan et al. 2013). However, recent genetic mapping studies using F 2 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 F 2 generation, and associations between transmission ratio distortion loci (TRDL) and quantitative traits (Chapman et al. 2013(Chapman et al. , 2016Brennan et al. 2014Brennan et al. , 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 F 2 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 ), S. squalidus in the UK is a stabilized hybrid descendent that has adapted to a new environment ). 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.

Mapping families
F 2 mapping families were founded from each of three paired cross combinations between three glasshouse grown (F 0 ) individuals representing each of S. aethnensis, S. chrysanthemifolius, and S. squalidus. Senecio aethnensis and S. chrysanthemifolius F 0 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 F 0 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 F 1 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 F 1 individuals with distinct maternal cytoplasms (i.e. F 1 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 F 1 ' s was chosen to found each of three F 2 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 et al. (2014) 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 et al. (2009) from all F 0 , F 1 and F 2 plants. The number of F 2 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 et al. (2014), 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 F 2 mapping families were genotyped for a total of 75 codominant genetic markers that were found to be polymorphic in the F 0 and F 1 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 et al. 2008), and 14 other codominant SSRs and indels derived from previously published Senecio sequences as described in Brennan et al. (2014). 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 et al. (2009).

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 et al. 2014). Because each F 2 mapping family was derived from two outbred F 0 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 cen-tiMorgans (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 et al. 1991, Fishman et al. 2001, map coverage according to Fishman et al. (2001), and genetic marker clustering using dispersion tests (Brennan et al. 2014).

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 100fold 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 F 0 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 F 0 parent, deficits or excesses of heterozygotes in terms of F 0 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 et al. (2016). Briefly, the genetic map was divided into 'n' intervals of equal length and the frequency of TRDL occurrence and cooccurrence 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 maplevel 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.

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, F 0 genotypes and genetic mapping information are provided in Supporting Information- Table S1. The full F2AS, F2AC and F2CS genetic maps are illustrated in Supporting Information- Figure S1 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 et al. 2014 and  LGs. See Supporting Information- Figure S1 for a depiction of all LGs corresponding to 10 chromosomes.  Linkage groups are named as in Brennan et al. (2014) 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). 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.

Genomic rearrangements
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 Supporting Information- Figure S1 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.

Transmission ratio distortion
Results of all TRD tests for genetic markers analysed in the F2AS and F2CS families are presented in Supporting Information-Tables S2 and S3, respectively. Equivalent results for the F2AC family are available in Table S2 of Brennan et al. (2014). The three F 2 mapping families differed significantly in the frequency of genetic markers showing genotypic TRD [Test 1, see Supporting Information- Table S4], 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 Supporting Information- Table S4]. 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 et al. 2014). 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.
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), # 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.
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 et al. 2014).

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 Supporting Information- Figure S2]. 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 Supporting Information- Figure S2]. No evidence was found for TRDLs being co-located with the midpoints of genomic rearrangments as inclusion of these data into the colocation analyses did not change the significance of test statistics [see Supporting Information- Figure S2].

Discussion
The genetic mapping and TRD analyses conducted on F 2 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 et al. 2015). 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

Genetic maps
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 et al. 2014). This finding was further supported by the results of an independent genetic mapping study conducted by Chapman et al. (2016). 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 Supporting Information-Fig. S1, 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 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). 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 et al. 2005;Ungerer et al. 2006) and Aegilops (Senerchia et al. 2016) 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 et al. 2013).
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 et al. 2011). 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 et al. 2004) and in Iris (Tang et al. 2010;Taylor et al. 2013). 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 et al. 2013;Osborne et al. 2013). 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 et al. 2010;Taylor et al. 2013), 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 et al. 2000;Rieseberg 2007, Karrenberg et al. 2007).

Transmission ratio distortion
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 et al. 2001;Harushima et al. 2001;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 Supporting Information- Table S4). 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 et al. 2013;Brennan et al. 2014;Chapman et al. 2016;Filatov et al. 2016), 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 F 2 mapping families involved full-sib F 1 crosses raising the possibility that some of the observed TRD could be due to biparental inbreeding. These species are self-incompatible and S. squalidus has been shown to suffer from inbreeding depression when selfed (Brennan et al. 2005(Brennan et al. , 2013. Inbreeding depression would affect patterns of TRD in the form of selection against reconstituted homozygous F 0 parental genotypes in the F 2 progeny. However, at TRDLs where tests of TRD of heterozygosity of marker alleles of different F 0 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 F 2 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 (Supporting Information- Table S4). 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 et al. 2005;Scascitelli et al. 2010;Senerchia et al. 2016;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 et al. 1995;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 et al. 2012;Filatov et al. 2016). 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 et al. 2014), 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 et al. 2007;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 et al. , 2014Chapman et al. 2016), 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 F 2 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 et al. 2014;Schumer et al. 2015b).

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 et al. 2017). 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 et al. 2017). 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.

Supporting Information
The following additional information is available in the online version of this article: Table S1. Summary of codominant molecular markers used for mapping. Table S2. Frequencies and segregation tests of genetic markers genotyped in the F2AS mapping family. Table S3. Frequencies and segregation tests of genetic markers genotyped in the F2CS mapping family. Table S4. Summary of genetic markers showing TRD on the F2AC, F2AS and F2CS genetic maps. Table S5. Summary of tests for co-location of TRDLs and rearrangements for each paired genetic map comparison. Figure S1. Interleaved genetic maps of all linkage groups from the F2AC, F2AS and F2CS mapping families. Figure S2. Probability of co-location of TRD loci and rearranged genomic regions for each paired genetic map comparison across a range of genomic size intervals.