From Hybrid Swarms to Swarms of Hybrids

Science has shown that the introgression or hybridization of modern humans (Homo sapiens) with Neanderthals up to 40,000 YBP may have led to the swarm of modern humans on earth. However, there is little doubt that modern trade and transportation in support of the humans has continued to introduce additional species, genotypes, and hybrids to every country on the globe. We assessed the utility of species distributions modeling of genotypes to assess the risk of current and future invaders. We evaluated 93 locations of the genus Tamarix for which genetic data were available. Maxent models of habitat suitability showed that the hybrid, T. ramosissima x T. chinensis, was slightly greater than the parent taxa (AUCs > 0.83). General linear models of Africanized honey bees, a hybrid cross of Tanzanian Apis mellifera scutellata and a variety of European honey bee including A. m. ligustica, showed that the Africanized bees (AUC = 0.81) may be displacing European honey bees (AUC > 0.76) over large areas of the southwestern U.S. More important, Maxent modeling of sub-populations (A1 and A26 mitotypes based on mDNA) could be accurately modeled (AUC > 0.9), and they responded differently to environmental drivers. This suggests that rapid evolutionary change may be underway in the Africanized bees, allowing the bees to spread into new areas and extending their total range. Protecting native species and ecosystems may benefit from risk maps of harmful invasive species, hybrids, and genotypes.


Introduction
By way of definition, the human species may be an instructive example of a hybrid swarm: the interbreeding (or introgression) of Homo sapiens and Neanderthal (Homo neanderthalensis) and other archaic human forms may have led to gains in disease immunity and cold-hardiness that helped modern humans spread to all corners of the globe [1]. Recent DNA analysis confirms small, but significant introgression of these humanoids in Asia, Europe, and Central Africa as late as 30,000 years ago [1][2][3]. In terms of species distribution modeling [4,5], it seems reasonable to assume that the inclusion of Neanderthal genes for cold-hardiness or disease resistance may have helped expand the "potentially suitable habitat" for modern humans. The subsequent swarms of modern humans, in turn, may have carried many other species in their travels [6][7][8].
Without a doubt, the increasing mobility of modern humans greatly diminished the former geographic isolation of other species. Modern trade and transportation, which has increased 40-fold in the past half-century, has significantly eliminated the geographic barriers that once separated the continents and the taxa on them [6,9]. As a result, many harmful, invasive plants, animals, and diseases are being spread across the globe with alarming consequences [8]. We focus on two areas of study receiving limited attention including: (1) the direct invasion of hybrid species, such as tamarisk (saltcedar; Tamarix) and Africanized honey bees, a genetic mixing of African and European honey bees; and (2) the potential for newly establishing alien species to hybridize with native or earlier-arriving alien species resulting in a hybrid swarm (i.e., a rapidly spreading hybrid taxon).
There are five urgent reasons to focus attention on hybrids. First, it has long been recognized that hybridization between native and introduced wildlife species already is a serious threat to conservation [10], with several cases of hybrids found in ducks, deer, fish, songbirds, and other species of management concern. Metcalf et al. [11] showed that "greenback cutthroat trout, Colorado River cutthroat trout, and rainbow trout readily hybridize, resulting in a hybrid swarm and the loss of the integrity of the cutthroat trout gene pool." In fact, historic fish stocking may have eliminated pure strains of native Colorado greenback cutthroat trout (Oncorhynchus clarki stomias) from their historic range [12]. Secondly, evidence of hybrid swarms is mounting. A hybrid of eastern and western Spartina, a cordgrass, spread extremely rapidly in the San Francisco Bay, negatively affecting wildlife habitat [13]. The hybrid females of 312 From Hybrid Swarms to Swarms of Hybrids domestic pigs and the introduced European wild boar are longer in length than either parent, and it is this hybrid that is spreading rapidly in the United States [14]. Thirdly, hybridization is more common in plants and animals than previously thought [15]. At least 25% of plants hybridize; and 10% of animals hybridize [16]. Fourthly, invasion can trigger rapid hybrid speciation [17], as evidenced by two native fruit flies hybridizing (Rhagoletis mendax x Rhagoletis zephyria ) in the northeastern United States to take advantage of an exotic nectar source, Japanese honeysuckle (Lonicera japonica) [18]. Lastly, with trade and transportation introducing taxa around the world at alarming rates, we agree with Vellend et al. [19] that ecologists urgently need to start mapping the distributions of new lineages.
Our goals are to: (1) provide preliminary maps and models of example hybrid taxa (or introgressing hybrid swarms); and (2) show the utility of commonly used species distribution models to raise awareness of the potential for hybrid swarms to spread. Common sense tells us we should document which taxa are likely to hybridize, produce copious fertile offspring, can disperse broadly, may avoid enemies better, or may cause significant economic or environmental harm or threats to human health. We think it's important to determine which hybrids may become a habitat generalist with a broader potential geographic range than their parents. This is not the definitive species distribution modeling paper. We acknowledge the shortcomings of small sample sizes, sample bias, and presence-only models [5,20,21], our imposed logistical defaults given the expense of genetics field and laboratory work. Instead, we present preliminary model results as "mapped hypotheses" to be improved with additional field work and modeling [22].
We draw on two examples where limited genetics information and occurrence data were readily available: Tamarix (a riparian shrub/tree) [23]; and Africanized honey bees [24]. Our specific objective was to demonstrate the utility of species distribution models to identify a potential hybrid swarm, such that the hybrid contributes significantly to the invasion process. This might be demonstrated if: (1) the hybrid is found to have broader potential habitat (i.e., more suitable habitat) than the parent taxa; or (2) if different genotypes of a species respond differently to environmental drivers leading to a greater extent of invasion.

Tamarix (saltcedars)
Saltcedars (Tamarix ramosissima Ledeb. and T. chinensis Lour.; family Tamaricaceae) are diploid deciduous shrubs or small trees. The species are distinct in their home lands: T. ramosissima is widely distributed across temperate Asia, while T. chinensis is restricted to China, Korea, and Japan [25]. Both species were introduced into North America in the 1800s, primarily in riparian habitats [26]. The genus is now widely spread across the western United States [27,28]. Recent multilocus DNA analysis of Tamarix by [29] found a much higher incidence of hybridization (83-87%) than previously reported and precipitated our current research. In this preliminary modeling project, we used locations of each taxon (T. ramosissima n=32 , T. chinensis n=29, and T. ramosissima X T. chinensis, n=32) to provide an indication of a hybrid swarm early in the invasion process.
We modeled the three Tamarix taxa with Maxent [4], a non-linear species distribution model that relies on presence data and drawing associations between what environmental conditions are available versus the environment where the species occurs. Preliminary Maxent models were run using the Software for Assisted Habitat Modeling version 1.1 [30], a tool that maintains a record of modeling parameters used and associated pre-processing steps to expedite and document the habitat modeling process. In this case, and in the case of many early invaders, it is difficult to differentiate absence due to unsuitability of habitat or absence due to a lack of propagules reaching the site, so a model that does not use absence data is the preferred alternative [31,32]. In Maxent, we set the maximum iterations to 5,000 to ensure convergence, ran 20 replicates withholding a subsampled 10% of the training data to test each replicate, and used a regularization value of 2 to control over-fitting. All other Maxent settings were left at defaults. We followed the target background approach, using presence of targeted genotypes as absence to limit sampling bias. Due to the low number of occurrence points and high correlation coefficients between pairs of variables, we restricted the number of environmental variables to two; minimum temperature of the coldest month and annual precipitation are known to be important predictors of Tamarix distribution at this scale [28]. This would also ensure 10 responses per independent variable [33]. Species distribution models, such as Maxent, fit correlations between the sampled distribution of an organism and environmental covariates such as climate factors or surrogates for primary productivity. To construct binary classifications of model output, we used the maximum of sensitivity plus specificity threshold. We used the latest version of freely available Maxent software, version 3.3.3k [4].

Apis (Africanized honey bee)
Apis mellifera includes approximately two dozen subspecies worldwide that are further classified into distinct lineages [34]. We classified A. mellifera subspecies into Africanized honey bees and European honey bees. We felt justified in this classification scheme based on a combination of literature review, previous swarm incidents reported in the media, and personal communication with bee biologists [24]. We further classified the Africanized honey bees into two distinct genetic groups based on mitotypes. Africanized honey bees are a genetic hybrid cross of Tanzanian A. mellifera scutellata and a variety of European honey bee including A. m. ligustica [35] that have been spreading north Environment and Ecology Research 2(8): 311-318, 2014 313 in the Americas since their introduction to Brazil in 1957. These hybrids first reached the United States from Mexico in 1990 and have continued their northward spread across the Southeastern United States [36]. Two suites of species distribution models were run with data on Apis mellifera lineages. We ran four different statistical modeling techniques: generalized linear models (GLM), boosted regression tree, multivariate adaptive regression splines, and random forest in SAHM [30] for Africanized bees (n = 573 presence; n = 264 absence points) and European honey bees (n = 253 presence; n = 584 absence points). These models were chosen for bee lineages because the presence or absence of dominant Africanized queens can be used in presence-absence models of Africanized and European bees. In other words, if an Africanized honey bee queen is present, then it may be assumed that a European honey bee queen is absent, and vice-versa. Following the methods in Jarnevich et al. [24], we initially considered 40 environmental data layers consisting of climate, land cover, and vegetation phenology variables to parameterize the models. Climate data included 19 bioclimatic layers from WorldClim [37], that are derived by interpolation of average monthly climate data at 30 arc second resolution. We included vegetation cover from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Continuous Fields (VCF) product, including percent estimates for trees (for 2005), herbaceous vegetation, and bare ground cover (for 2001) [38]. The MODIS Land Surface Phenology product provided 15 metrics of seasonal variation in vegetation productivity from 2001-2007, excluding 2005 [39]. Similar to the long term average climate used, we calculated the average of each phenology metric across all years available at the time of analysis (2001-2007 excluding 2005). All data were resampled to 30 arc seconds (~1km) to match the finest spatial resolution data set in ArcGIS v. 10.1 as detailed in Jarnevich et al. [24].
We tested for cross-correlation among the 40 environmental data layers in SAHM using the maximum of Pearson, Spearman, and Kendall correlation coefficients. When two or more variables were highly correlated (r ≥ 0.70 or ≤ -0.70), variables with lower biologic relevance were dropped (following Jarnevich et al. [24]). Each of the four statistical models began with the same set of environmental variables (Table 1); however, each algorithm except random forest incorporated variable selection, therefore the final set of variables used in each model was slightly different. The GLM employed standard stepwise regression using Akaike's information criterion. We used 10 fold cross-validation to assess model performance.
After all four statistical models were run for European and Africanized honey bees, we combined the SAHM [30] continuous probability surface outputs in ArcGIS v. 10.1 to produce ensemble maps of the predictions. The sensitivity equals specificity threshold for each model was used to classify each grid cell as suitable or unsuitable [40]. This threshold rule ensures presence locations are equally likely to be false positives as absence locations. The resulting binary maps for each model were then added, resulting in a map that displayed areas where one, two, three, or all four models agree [24]. Finally, to evaluate the "tension zone", where habitat suitability for European and Africanized honey bees currently overlap, we reclassified the probability surface for each model into a binary map using (threshold value +/-0.10), combined the four models into an ensemble for each lineage, and then further combined the two ensemble maps into a new ensemble where values of zero, one, or two (number of models in agreement, respectively) equal no tension and values of three or four equal tension where at least three of the statistical models predict European and Africanized honey bees overlap in habitat suitability.
Second, mitochondrial DNA analysis on 89 Africanized bees allowed us to model the potential distribution of 72 occurrences of the A1 mitotype (common in Africa, Brazil, and Mexico), and 17 occurrences of the A26 mitotype (common in SW Africa, but rare in the Americas), using the same Maxent approach that we used for Tamarix modeling above, and, to avoid over-fitting the model, we again selected only two variables (minimum temperature of coldest quarter and season length) as predictors. We used other AHB locations as background following the target background approach [41]. These environmental variables have been shown to be important drivers of Africanized bee distribution [24]. Similar to models described above, we used the maximum of sensitivity plus specificity as the threshold for binary classifications. In all our models, we recognized that our results may have been affected by small sample size, but Maxent has worked well with small sample sizes [42,43]. Our results also may be affected by sample location bias (a point we return to in the discussion).
Regardless of the species, we used the area under the receiver-operating characteristic (ROC) curve (AUC) to assess model accuracy with each validation dataset. AUC, a measure of model fit, has been widely used in several model comparison studies. Usually AUC values of >0.9 indicate high accuracy, values of 0.7-0.9 indicate good accuracy, and values 0.5 (random) to 0.7, indicate low accuracy. We minimized the potential problems associated with the reliance on AUC values [44] by: (1) selecting models appropriate for the data (e.g., presence or presence-absence models); (2) selecting only models that performed well in previous model comparison studies [45]; and (3) eliminating one of each pair of highly cross-correlated predictor variables prior to modeling (i.e., in the honey bee models). We also used percent correctly classified (threshold dependent metric where sensitivity equals specificity) as an indication of model performance.

Tamarix
Preliminary Maxent models were reasonably strong (test AUC > 0.83) for Tamarx ramosissima, T. chinensis, and the hybrid T. ramosissima x T. chinensis, despite the modest sample size of about 30 observations per taxa (Fig. 1). We were surprised that the extent of the potential habitat suitability for the hybrid T. ramosissima x T. chinensis hybrid was slightly greater than the parent taxa (T. ramosissima or T. chinensis; Fig. 1). Selected response curves from the Maxent models for minimum temperature of the coldest month and annual precipitation showed that the hybrid T. ramosissima x T. chinensis response curves were intermediate between T. ramosissima and T. chinensis (Fig. 2), describing a more general response to the environmental variables.

Africanized Honey Bees
All statistical models produced test AUC values greater than 0.85 and correctly classified suitable habitat greater than 74%, indicating strong performance ( Table 1). The models show that Africanized honey bee and European honey bees have diametrically opposite areas of habitat suitability and occurrence (Fig. 3). We overlaid binary probability maps of the Africanized honey bee and European honey bee distributions to examine the zones of tension in the two taxa (Fig. 4). Each map was calculated by reclassifying the continuous model prediction to be ±10% around the binary threshold value for each model, which identifies the areas of moderate suitability. The western Great Basin, western Great Plains, and southern Appalachian mountains are regions highlighted by these models, where the European honey bees may be struggling to maintain dominance. Maxent models of the two dominant mitotypes of the Africanized honey bees were the strongest models of this study (AUC ~0.90; Fig. 5). Selected response curves for minimum temperature of coldest month and annual precipitation suggest that the two mitotypes may respond differently to environmental drivers. The A26 mitotype is associated with slightly cooler and drier sites compared to the more common A1 mitotype (Fig. 5)

Caveats
We realize that we have small datasets for the Tamarix genotypes, and for the mitotypes of the Africanized honey bees, due largely to the high costs of obtaining genetic data. Small sample size and sample location bias may have affected our model results, but sample sizes as low as n=10 with Maxent performed well [42,43]. Our purpose here is to demonstrate the utility of modeling hybrid swarms and wildly introgressing populations to provide preliminary risk maps as "mapped hypotheses" to be improved with additional field work and modeling [22]. Since the various genotypes of Tamarix and mitotypes of bees were modeled in similar ways, we think they serve our purpose..

316
From Hybrid Swarms to Swarms of Hybrids

Potential for Hybrid Range Complementarity
Our preliminary investigation of Tamarix hybrids and Africanized honey bees yielded a similar, frightening result: at least some hybrids may be highly invasive by: (1) responding differently to environmental drivers than the parent taxa (Figs. 2 and 5); or (2) occupying slightly different habitats than the parent taxa ( Figs. 1 and 5). This complementarity of taxa distributions may contribute to the spread, extent, and persistence of hybrid taxa. Other mechanisms may also be important for hybrid persistence such as avoiding pathogens and other adaptations [46].

Conservation Concerns of Hybrid Swarms
We provided two recent examples of hybrid swarms that are not without conservation concerns. Friedman et al. [47] and Gaskin and Kazmer [29] found that the F1, F2, and backcrosses to two parent species, T. ramosissima and T. chinensis, may represent 83-87% of the tamarisk invasion in the western United States. Recent studies have found that this complicates the biological control of the genus, because increased levels of T. ramosissima introgression resulted in higher investment in roots and tolerance to defoliation, and less resistance to biological control agents often used by managers to decrease the range and spread of this genus [48]. Plant hybridization, particularly in this harmful invasive species, may hinder containment efforts in natural areas.
Likewise, the hybridized Africanized honeybees may have many evolutionary advantages over European honeybees and native bees. The Africanized honeybees have faster growth rates than European honeybees, resulting in more swarms to dominate an area [49]. And, since European honeybee queen mate disproportionately with African drones, this is likely to cause rapid displacement of European honeybee genes in a colony. For example, when a queen is inseminated with an equal portion of African drone semen and European honeybee semen, "the queens preferentially use the African semen first to produce the next generation of workers and drones, sometimes at a ratio as high as 90 to 10" [49]. According to bee expert, David Roubik [50], Africanized honeybees "groom themselves more often than Italian bees, making them less likely to get sick from mites and other parasites" that plague European honeybees. In South America, the Africanized honeybees withstand broader environmental gradients, from desert to rainforest. In the more extreme environments "the European bee will just starve to death, the Africanized bee is foraging, bringing things in, and getting by," said Roubik. Combined with their aggressive behavior, these traits and characteristics may explain why Africanized honeybees are displacing native bees in Mexico and elsewhere [51].
We in Florida [52]. Vilà and D'antonio [53] demonstrated hybrid vigor for clonal growth in an invasive alien Carpobrotus (Aizoaceae) in coastal California. The native California Tiger Salamanders (Ambystoma californiense) hybridized with the introduced Barred Tiger Salamanders (Ambystoma tigrinum mavortium) resulting in a more robust hybrid [54].
As global trade increases, we expect many more closely related species to come in contact [55], setting the stage for swarms of hybrids and rapid evolution. We echo the concerns of Lee [56] more than a decade ago that conservationists should "emphasize the utility of genomic approaches for determining invasion mechanisms, through analysis of gene expression, gene interactions, and genomic rearrangements that are associated with invasion events. "The invasion potential of better adapted hybrid swarms may increase with increasing levels of trade and transportation allowing for congeners from other countries and regions to come within mating distance for the first time in a long time, or at any time in the past [16,17]. Introgression may play an increasingly important role in the distribution of taxa [13,16]. For example, a cline of hybridization, cold-hardiness, and introgression of taxa in the genus Tamarix has been observed along a latitudinal gradient in the United States [47], suggesting rapid adaptation may be possible and expanding the potential invaded range. Likewise, multiple introductions of Africanized honey bees add a complexity in analyzing invasion patterns [24]. Our findings indicate that honey bee lineages differ in their climate associations and genotypic niche modeling may offer insight into the invasion dynamics of these hybrids.
It is highly likely that trade and transportation, in support of our hybrid-human swarm, will continue to introduce additional species, genotypes, and hybrids to every country on the globe [8]. If our conservation goal remains to protect native species and ecosystems [57], risk maps and distribution models of harmful invasive species, hybrids, and genotypes may help [58,59].