Inferring responses to climate warming from latitudinal pattern of clonal hybridization

Abstract Climate warming may affect reproductive isolation between sympatric sister species by modifying reproductive phenology or mate choice. This is expected to result in a latitudinal progression of hybridization in response to the shifting of environmental conditions. The fish species northern redbelly dace (Chrosomus eos) and finescale dace (C. neogaeus) display a wide sympatric distribution in North America. The asexual reproduction of their hybrids allows determining where and when hybridization occurred. The aim of this study was twofold: first, to assess whether temperature affected reproductive isolation, and second, whether the effects of climate warming resulted in a latitudinal progression of hybridization. We performed a 500 km latitudinal survey (51 sites) in southeastern Quebec (Canada) and determined the distribution of clonal hybrid lineages. Results revealed a total of 78 hybrid lineages, including 70 which originated locally. We detected a significant difference between the southern and northern range of the survey in terms of the proportion of sites harboring local hybrids (20/23 vs. 8/28 sites, respectively) and hybrid diversity (57 vs. 13 lineages, respectively). This confirmed that there was more frequent interspecific mating in the warmest sites. In the southern range, diversity of lineages and simulations suggest that hybridization first took place (>7,000 years) in sites characterized by a longer growing season, followed by northerly adjacent sites (ca. 3,500–5,000 years). Moreover, evidence of hybridization occurring in present‐day time was detected. This suggests that the current warming episode is going beyond the limits of the previous warmest period of the Holocene.


| INTRODUC TI ON
Global warming refers to the general increase in temperature across the planet since the late 19th century. The modifications of environmental conditions brought by the rise of the temperature have major impacts on premating barriers, making hybridization a serious consequence of global warming (Chunco, 2014;Hoffmann & Sgro, 2011;Parmesan & Yohe, 2003;Scheffers et al., 2016). Firstly, shift in the geographic distribution of species is one response to the current climate warming (Chen, Hill, Ohlemüller, Roy, & Thomas, 2011;Kelly & Goulden, 2008) and also occurred during postglacial expansion (e.g., Dyke, 2005). Secondary contact resulting of a differential shift in the range of closely related species may, however, lead to the formation of a hybrid zone and have been reported for several examples (e.g., Ryan et al., 2018). Secondly, climate warming can also affect reproductive phenology or mate choice and weaken reproductive isolation between sympatric nascent species (reviewed in Chunco, 2014). The reproductive phenology of organisms is strongly influenced by temperature, and warmer conditions often result in earlier time of mating (Benard, 2015;Carey, 2009;Walther et al., 2002;Wegge & Rolstad, 2017). Differential shift in the timing of reproduction may lead to an overlapping of reproductive periods between species and to hybridization (Chunco, 2014;Vallejo-Marín & Hiscock, 2016). Warming can affect mate choice when habitats become unfavorable to one species.
This may increase the probability of hybridization due to asymmetric availability of mates. Indeed, according to the Hubbs principle (Hubbs, 1955), a species decreasing in abundance is expected to breed with individuals of another more abundant species as conspecific individuals become scarce (Avise & Saunders, 1984;Grant & Grant, 1997;Jansson, Thulin, & Pehrson, 2007;Randler, 2002;Wirtz, 1999).
The effect of warming is not uniform across the latitudinal range of species because temperature varies with latitude (and altitude). If temperature must reach a threshold to trigger a given change, warmer regions will be first affected, and with time (and continuous warming), the geographic limit reaching this threshold will move latitudinally. The modification of premating barriers between sympatric related species is therefore expected to result in a wave of hybridization through the sympatric area in response to climate warming.
Moreover, the current climatic warming is not an isolated event.
The transition from Pleistocene to Holocene ca. 12,000 years BP was marked by a strong increase in temperature, while the following Kyears were punctuated by cooling and warming episodes (Marcott, Shakun, Clark, & Mix, 2013). For instance, the Holocene Climate Optimum was a period characterized by the warmest temperatures since the end of the Pleistocene, with conditions warmer than now (but see Porter et al., 2019). The time of Holocene Climate Optimum varied widely across the planet (Kaufman et al., 2016;Marcott et al., 2013). For instance, this peak warming occurred ca. 6,000-3,500 years BP in northeastern North America (Comtois, 1982;Larochelle, Lavoie, Grondin, & Couillard, 2018).
To determine the extent of historical warming events and assess the magnitude of the current episode on hybridization, the clonal hybrids represent a useful biological indicator through which to infer the effects of climate warming since they keep track of hybridization events.
The northern redbelly dace (Chrosomus eos; Figure 1) and the finescale dace (C. neogaeus; Figure 1) are fish species with a wide overlapping distribution in the north of North America, mostly located in regions previously covered by the last Pleistocene ice sheet (Scott & Crossman, 1973). Distinct spawning times are expected to provide a premating barrier between those species (Das & Nelson, 1990). However, this temporal barrier is not completely impervious and C. eos-neogaeus hybrids are also widely distributed (Angers & Schlosser, 2007). A particular characteristic of these hybrids is that they are only females and reproduce clonally, as observed in several other unisexual hybrid vertebrates (Neaves & Baumann, 2011). The absence of meiosis in the hybrid genome prevents backcrosses and genetic exchanges between C. eos and C. neogaeus species (Lafond, Hénault, Leung, & Angers, 2018). Therefore, once hybrid individuals are produced, their genotypes can persist through time via clonal lineages.
Surprisingly, most of those hybrids do not originate from local hybridization, but rather predated the end of the Pleistocene.
Indeed, previous studies have reported the presence of a single or a few widespread lineages in most of the regions surveyed (Angers & Schlosser, 2007;Goddard, Dawley, & Dowling, 1989). However, both species can still interbreed and local hybridization events have been genetically confirmed in the southeastern region of Quebec, Canada (Vergilino, Leung, & Angers, 2016). Compared to Pleistocene migrants, local hybridization resulted in several lineages specific to a given site (private lineages). These assumptions parallel those inferred from mitochondrial DNA in phylogeographic studies of regions formerly covered by glaciers (e.g., Bernatchez & Wilson, 1998). mtDNA is asexually transmitted as are the clonal hybrid lineages. On the one hand, mtDNA haplotypes already present in glacial refuges displayed a large geographic distribution. They took advantage from F I G U R E 1 The northern redbelly dace (Chrosomus eos) and the finescale dace (C. neogaeus). Clonal hybrids C. eos-neogaeus (absent in the picture) result from hybridization events between those species. The northern redbelly dace has two black stripes alongside (right side of the picture), while the finescale dace exhibit a thick orange stripe (left side of the picture) the proglacial lakes and temporary bridges among hydrographic networks to spread throughout different watersheds (Curry, 2007;Mandrak & Crossman, 1992;Parent & Occhietti, 1999). On the other hand, mtDNA variants that appeared more recently are generally restricted to the contemporary drainage system as interconnections are no longer present.
Since local hybridization events appear to be rare in most sperm-dependent unisexual hybrids, this raises questions about the causes of such high rates of spontaneous hybridization in that region. The aim of this study was twofold: first, to assess whether temperature (or other correlated environmental conditions) affected reproductive isolation between C. eos and C. neogaeus (spatial pattern), and second, whether effects of climate warming (current and/ or historical) resulted in a latitudinal progression of hybridization through time (temporal pattern). To address these objectives, we performed a 500 km latitudinal survey in the southeastern Quebec region to investigate the large-scale geographic variation in distributions of the C. eos-neogaeus hybrids, as a consequence of local hybridization between C. eos and C. neogaeus.
Clonal hybrids are not only useful to assess local hybridization that occurred thousands of years ago, but they could also be used to infer time since hybridization. This would be possible if hybridization is not a continuous process, so in the absence of new hybrids, lineage diversity decreases with time. Complete extirpation of one parental species stops hybridization. Such a scenario is one conclusion of a species decline associated with hybridization according to the Hubbs principle (Hubbs, 1955). However, the limited occurrence of local hybridization when both species co-occur suggests that other factors may prevent hybridization.
Chrosomus neogaeus first reproduces in early spring, while C. eos reproduces later (Das & Nelson, 1990). The hybrid lineages analyzed so far were always the result of reproduction between late C. neogaeus females and early C. eos males, since all hybrids detected so far harbor the maternally inherited mitochondrial DNA of C. neogaeus (Angers & Schlosser, 2007;Binet & Angers, 2005;Goddard et al., 1989;Vergilino et al., 2016). Hybrids are known to display an intermediate reproduction period between that of C. neogaeus and C. eos. This reproduction period must overlap with that of one of the parental species because the reproduction of the hybrids requires the sperm of a parental species to activate the development of their eggs (gynogenesis). In addition, the asexual reproduction of all-female clonal hybrids is expected to result in a rapid demographic expansion of hybrid population . The high abundance of sperm-dependent female hybrids compared to the low number of remaining late C. neogaeus females may prevent or strongly reduce the probability of further hybridization events.
If hybridization is a punctual event rather than a continuous process, the number of hybrid lineages is expected to decrease with time in the absence of new hybridization events. A single interspecific mating generates a high diversity of lineages since each zygote producing a female can evolve as a distinct hybrid lineage (no androgenesis was reported in this species complex). However, through generations, random processes related to lineage sorting, as well as natural selection, progressively reduce the number of lineages. Therefore, the diversity of lineages at a given site depends on the time elapsed since the hybridization events (in addition to population size and environmental conditions). Relative Note: List of sampled sites with their geographic coordinates and sample size of each biotype (E, N, and H refer to Chrosomus eos, C. neogaeus, and the hybrids, respectively). The number of private hybrid lineages (and shared lineages) of each site is indicated. Lineages shared among sites are indicated by the same letter.
LGS refers to the length of the growing season in days according to Agriculture and Agri-Food Canada (2014). Sampling refers to sites sampled 1- Vergilino et al. (2016) and 2-this study. Each site is designated by two letters referring to its hydrographic networks. Sites with more than one number (e.g., RI 2-4) indicated sites pooled together.

TA B L E 1 (Continued)
time since local hybridization events can be correlated to lineage diversity.
In this contextual framework, we can formulate specific predictions for each objective. First, we expect a higher occurrence of hybridization in southern sites than in northern sites due to the higher temperature at lower latitude (closer to the threshold expected to weaken premating barriers). In northern sites, where premating barriers are not yet disturbed, the presence of locally produced hybrids is expected to be infrequent, even in the presence of sympatric parental species. Second, if climate warming (previous episodes during Holocene as well as the current episode) modifies premating barriers, hybridization is expected to progress northward through time.
We therefore predict that the hybrid diversity resulting from local hybridization will vary according to a latitudinal gradient. Since the number of lineages decreases with time since hybridization, diversity is expected to increase going northward, as hybridization events are more and more recent.

| Sampling sites
We analyzed a total of 51 sites along a 500-km-long transect, from 45°02′ to 48°29′ north ( Figure 2, Table 1 integrative index of climate differences among regions (Linderholm, 2006) and is often correlated to phenotypic change in animals (Garel, Solberg, SAEther, Herfindal, & Høgda, 2006;Peñuelas, Filella, & Comas, 2002;Zani, 2008). Increasing temperature is expected to lengthen the growing season of temperature-limited organisms such as ectotherms. The LGS was measured in days, starting from 10 days after average daily temperature is above 5°C and ending when the minimum daily temperature is 0°C or October 31st is reached-whichever comes first and has been estimated from records, using a 29-year (from 1971 to 2000) reference period of measure (Agriculture & Agri-Food Canada, 2014). We used the map of Agriculture and Agri-Food Canada (2014) where LGS categories are incremented by 10 days, to assign our sampling sites to one or the other categories (Figure 2). A difference of 40 days in the LGS was observed between the southernmost and the northernmost sites.
Although the climatic conditions were likely different than those during the Holocene, these categories were used to emphasize the latitudinal differences along our 500 km survey.

| Diversity of hybrid lineages
External morphological characteristics were used to visually identify sampled individuals as C. eos, C. neogaeus, or hybrids (New, 1962).
DNA was extracted from a small piece of the upper lobe of the caudal fin, and we confirmed the visual identification using the genetic markers specific to nuclear and mitochondrial genomes developed by Binet and Angers (2005).
A hybrid lineage was defined by individuals originating from the fertilization of one C. neogaeus egg by a C. eos sperm. As the clonal reproduction of hybrids leads to the absence of recombination and segregation, individuals of a given lineage are expected to display the same multilocus genotype (Angers & Schlosser, 2007).
However, due to the high mutation rate of microsatellite loci (Estoup & Angers, 1998), variants could be detected. They generally differed by a single-step mutation at one or a few loci. A lineage was thereby characterized by a consensus genotype defined by the allele of the invariant loci and the most abundant allele for the variable loci (Angers & Schlosser, 2007).
In the Chrosomus eos-neogaeus complex, a high proportion of triploid hybrids are produced by the incorporation of the sperm genome into the diploid hybrids' eggs (Goddard et al., 1989;. Since triploid hybrids display both the consensus multilocus genotype of a given hybrid lineage and the spermatozoid haplome, it was possible to identify the hybrid lineage from which they derived. However, when no match could be made between a triploid genotype and any of the consensus lineages, this triploid genotype was considered as an additional lineage for which it was not possible to determine the hybrid genotype.

| Local hybridization versus pleistocene migrants
Once lineages were identified, the next step was to determine whether hybrids were migrants from Pleistocene populations or descendants of local hybridization. The occurrence of the same lineage through different hydrographic networks was considered as the signature of hybrid lineage from the Pleistocene (Angers & Schlosser, 2007;Vergilino et al., 2016). In contrast, hybrid lineages produced locally occurred in isolated hydrographic networks similar to those in the current landscape and were expected to display a narrow geographic distribution limited to distinct subdrainage basins (Vergilino et al., 2016). Lineages specific to a given site (private lineages) were therefore likely the result of local hybridization events.

| Geographic distribution of species and hybrids
The evenness of species distribution (in terms of the proportion of sites occupied) across LGS categories was tested by Fisher's exact test, while the abundance of species was compared by a chi-square test.

| Spatial pattern of hybridization
To assess the effect of temperature on reproductive isolation of C. eos and C. neogaeus, the difference in the number of sites harboring private hybrid lineages (as indication of local hybridization) and diversity of private lineages were tested according to the LGS using Fisher's exact test and nonparametric Mann-Whitney U tests, respectively. The diversity of private lineages across a latitudinal gradient was also represented using a LOESS/LOWESS (Locally Weighted Scatter-plot Smoother) regression (Cleveland, 1979;Cleveland, Grosse, & Shyu, 1991). The smoothing parameter was selected automatically using a generalized cross-validation criterion (Golub, Heath, & Wahba, 1979), as implemented for fANCOVA package version 0.5-1 (Wang & Wang, 2010) in the statistical programming environment R 3.4.0.

| Temporal pattern of hybridization
To assess whether hybridization occurred first in the southernmost sites and then progressed northwardly, we compared the diversity

| Geographic distribution of species and hybrids
Both C. eos and C. neogaeus, as well as hybrids, were found throughout all hydrographic networks ( Figure 3, Table 1). The sexual species C. eos was detected in most of the sites (40 out of 51 sites; 78.4%).
This species displayed an even distribution; the proportion of sites occupied by C. eos did not differ among LGS categories (Fisher's exact test p = .3745).
The other sexual species C. neogaeus was only detected in 18 sites (35.3%) and was less abundant than C. eos (χ 2 = 17.63; df = 1; p = .00003). The distribution of this species is uneven (Fisher's exact test p = .0429), being less abundant in the LGS of 150-160 days ( Figure 3). For instance, this species occupied a single site in the southernmost networks (YA, RI an SF; Table 1).
At the opposite of C. neogaeus, hybrids were less abundant in the northern sites (LGS of 120-130 days; Figure 3).

| Local hybridization versus pleistocene migrants
The genetic survey of hybrid individuals revealed a total of 78 distinct lineages (Table 1); 41 lineages were detected in the new sites surveyed in this study from CH, SJ, OU, LO, RM, and MI networks, in addition to the 37 lineages previously detected in the YA, RI, and SF networks (Vergilino et al., 2016). Although the sampling size differs among sites, there is no correlation between the number of hybrid individuals sampled and the number of lineages detected (R 2 = 0.08; p = .1).

| Spatial pattern of hybridization
Private hybrid lineages were present in 28 sites (54.9%). Their distribution across LGS categories is also uneven (Mann-Whitney

| Temporal pattern of hybridization
We assessed the temporal change, more specifically whether Simulations of the random processes related to lineage sorting clearly illustrated the reduction in the number of lineages through time, when hybridization is a punctual event rather than a continuous process ( Figure 5). Diversity remains high in the whole system after 100 generations (ca. 80 lineages). However, the number of lineages per population rapidly decreased to 50 lineages per population (the number expected in site CH-4) in less than 20 generations.
A system comparable to LGS of 140-150 days with a diversity of 5-6 lineages at some sites is detected when hybridization occurred between ca. 3,000 and 5,000 years BP. A system comparable to LGS of 150-160 days (the southernmost sites) with a diversity of three lineages or less is detected when hybridization occurred >7,000 years BP. Finally, hybridization that occurred 10,000 years ago (the maximum of generations corresponding to postglacial expansion) still displays diversity and some systems could maintain up to three lineages.

| D ISCUSS I ON
The aim of this study was first to assess whether temperature (or associated environmental conditions) affected reproductive isolation between C. eos and C. neogaeus (spatial pattern) and second whether climate warming (previous or current episodes) resulted in a latitudinal progression of hybridization through time (temporal pattern).
The clonal hybrids that keep track of hybridization between those species are therefore a useful biological indicator through which to infer the effects of climate warming.
Our 500 km latitudinal survey revealed the presence of 78 distinct C. eos-neogaeus hybrid lineages. This striking result contrasts with diversity from the other regions of North America surveyed so far. Previous studies reported a single or a very low number of widespread lineages per region (Angers & Schlosser, 2007;Doeringsfeld, Schlosser, Elder, & Evenson, 2004;Elder & Schlosser, 1995;Goddard et al., 1989;Vergilino et al., 2016). Such a low diversity in very large areas is typical of mtDNA diversity of species that colonized regions F I G U R E 4 Relationship between hybrid diversity and latitude. Scatter plot and local regression curve between the number of private lineages and latitude.

Lineages/site
LGS ( covered by the ice sheet during the Pleistocene (Bernatchez & Wilson, 1998). Applied to the context of clonal fishes, low diversity has been associated with postglacial expansion of a single (or a few) clones from glacial refuges (Angers & Schlosser, 2007).
Results of our study revealed that most of the hybrid lineages (70 lineages) were specific to a single site or were detected in a few geographically close sites that shared the same lineage assemblage.
The very high diversity of lineages and their restricted distribution are consistent with the hypothesis that hybridization events occurred locally. Genetic analyses provided another support to this hypothesis as the multilocus genotype of hybrids from RI, YA, and SF networks matched that of C. eos sympatric populations, one of the species involved in hybridization (Vergilino et al., 2016).  (Bernatchez & Wilson, 1998;Legendre & Legendre, 1984).
Altogether, these results support the rejection of the postglacial colonization and suggest that these 70 private lineages are the result from local hybridization. The next step is therefore to investigate local factors that affected the latitudinal variation in lineage diversity.

| Spatial pattern of hybridization
The spatial distribution of private hybrid lineages is consistent with a temperature or environmental condition effect on hybridization. While C. eos and C. neogaeus are (or were) sympatric in most of the hydrographic networks analyzed, the diversity of private hybrid lineages strongly varied latitudinally. We observed a higher proportion of sites where hybridization occurred locally and a higher number of private lineages in the southern half of the survey than in the northern sites.

| Effect of temperature on reproductive barriers
Hybridization may result from the decline of one species. When the abundance of conspecific mates decreases, the reproduction with more abundant heterospecific males is then favored (Hubbs, 1955).
This hypothesis finds support in the quasi-absence of C. neogaeus from the networks of southeastern Quebec (RI, YA, and SF;

| Temporal pattern of hybridization
We predicted that if local hybridizations were the consequence of climate warming, diversity of hybrids would vary latitudinally.
Diversity of lineages from local hybridization events is expected to decrease as the time elapsed since the hybridization events increase, in the absence of additional hybridization. The increase in the number of lineages detected when moving northward along the southern half of the survey is consistent with this prediction.
The high hybrid diversity in the LGS of 140-150 days reflects the most recent historical hybridization events observed. This region also harbored four sites in which at least five distinct private lineages were detected, a diversity of lineages observed in no other sites or previous studies (Angers & Schlosser, 2007;Doeringsfeld et al., 2004;Elder & Schlosser, 1995;Goddard et al., 1989). Simulations indicated that a system maintaining such diversity originated from hybridization that should have occurred between ca. 3,000 and 5,000 years ago. This scenario is consistent with the Holocene Optimum time expected to occur ca. 6,000 to 3,550 years BP in southern Quebec (Comtois, 1982;Larochelle et al., 2018).
In spite of more exhaustive sampling (in terms of a larger longitudinal survey) performed in the LGS of 150-160 days Vergilino et al., 2016), and a high number of lineages detected, the diversity of private lineages per site remains lower than in the northerly adjacent LGS category. This lower diversity suggests that hybridization took place earlier in the southernmost sites (RI, YA, and southern part of SF). Results of simulations indicated that hybridization in sites maintaining a maximum of three lineages likely occurred > 7,000 years BP. These events occurred long before (around 2,000 years) those observed in the sites of LGS of 140-150 days (northern part of SF and CH). While these simulations must be interpreted with caution, they provide a likely scenario illustrating the timing of hybridization events which occurred in southern Quebec during warming of the Holocene Climate Optimum period and the current warming episode.

| Effect of current global warming
One site (CH-4) revealed an exceptionally high number of lineages, one order of magnitude higher than in other sites. The unusually high diversity of CH-4 contrasts with that of the other sites of the survey. The low diversity of those sites is expected to be representative of the population. For example, the same lineages were detected when sampling geographically close sites (Vergilino et al., 2016), indicating that most of the lineages were recovered. While 15 hybrid lineages including 13 private lineages were sampled in CH-4, we estimated that this population should include at least 50 different lineages to produce diversity similar to our sampling. This suggests very recent (during the last decades) hybridization events occurred at CH-4 site.
Effects of the current climate warming are suspected because of a significant increasing trend in mean water temperature and precipitation, associated with increasing temperature, has been observed over the past decades in southern Quebec (Hudon, Armellin, Gagnon, & Patoine, 2010;Yagouti, Boulet, Vincent, Vescovi, & Mekis, 2008;Zhang, Vincent, Hogg, & Niitsoo, 2000). As evidence of warming, the mean annual temperature recorded between 1960 and 2005 increased by ca. 1.0°C in the LGS of 140-150 and 150-160 days, while no significant changes were detected in the northern part of the sampling (Yagouti et al., 2008).
Results of the simulations revealed that 100 lineages produced by a local hybridization can rapidly decrease to 50 lineages per population (the number expected in site CH-4) in less than 20 generations. This result suggests that a LGS of 140-150 days could serve as a reference for conditions triggering hybridization between the northern redbelly dace (C. eos) and the finescale dace (C. neogaeus).
In conclusion, the present study provides an empirical framework to assess how historical climate warming influenced hybridization between sympatric sister species. However, detection of contemporary hybridization close to the northern limit of historical hybridization suggests that the current warming episode is going This suggests that organisms will soon have to face effects of the current climate warming they never experienced in the last era.

ACK N OWLED G M ENTS
We are grateful to reviewers who provided thoughtful comments on earlier drafts of this paper, as well as to the numerous students who took part in this project. This research was supported by a research grant from the Natural Sciences and Engineering Research Council of Canada (NSERC) to BA (#238600).

DATA AVA I L A B I L I T Y S TAT E M E N T
All relevant data are within the paper.