Influence of invasion history on rapid morphological divergence across island populations of an exotic bird

Abstract There is increasing evidence that exotic populations may rapidly differentiate from those in their native range and that differences also arise among populations within the exotic range. Using morphological and DNA‐based analyses, we document the extent of trait divergence among native North American and exotic Hawaiian populations of northern cardinal (Cardinalis cardinalis). Furthermore, using a combination of historical records and DNA‐based analyses, we evaluate the role of founder effects in producing observed trait differences. We measured and compared key morphological traits across northern cardinal populations in the native and exotic ranges to assess whether trait divergence across the Hawaiian Islands, where this species was introduced between 1929 and 1931, reflected observed variation across native phylogeographic clades in its native North America. We used and added to prior phylogenetic analyses based on a mitochondrial locus to identify the most likely native source clade(s) for the Hawaiian cardinal populations. We then used Approximate Bayesian Computation (ABC) to evaluate the role of founder effects in producing the observed differences in body size and bill morphology across native and exotic populations. We found cardinal populations on the Hawaiian Islands had morphological traits that diverged substantially across islands and overlapped the trait space of all measured native North American clades. The phylogeographic analysis identified the eastern North American clade (C. cardinalis cardinalis) as the most likely and sole native source for all the Hawaiian cardinal populations. The ABC analyses supported written accounts of the cardinal's introduction that indicate the original 300 cardinals shipped to Hawaii were simultaneously and evenly released across Hawaii, Kauai, and Oahu. Populations on each island likely experienced bottlenecks followed by expansion, with cardinals from the island of Hawaii eventually colonizing Maui unaided. Overall, our results suggest that founder effects had limited impact on morphological trait divergence of exotic cardinal populations in the Hawaiian archipelago, which instead reflect postintroduction events.


| INTRODUC TI ON
The recognition that biological invasions provide unique insight into the mechanics of evolutionary divergence has led to a spike in published research on postestablishment evolution of exotic species (e.g., Dlugosch & Parker, 2008a,b;Suarez & Tsutsui, 2008). There are now several examples of marked divergence in genetic or phenotypic traits between two or more exotic populations (Egizi, Fefferman, & Fonseca, 2015;Freed, Conant, & Fleischer, 1987;Lucek, Sivasundar, & Seehausen, 2014;Phillimore et al., 2008;Westley, Conway, & Fleming, 2012;Xu et al., 2010). Such differences can be explained by in situ adaptation to local biological and environmental conditions, or from events that occurred within the species' invasion history (e.g., Allendorf & Lundquist, 2003;Dlugosch & Parker, 2008a;Keller & Taylor, 2008). In particular, founder effects can result in divergence of traits across exotic populations if colonizing individuals are derived from two or more genetically and/or phenotypically structured native subpopulations and introduced in such a way where these features are structured across the exotic range (Keller & Taylor, 2008). Here, we deduce the invasion history of northern cardinals (Cardinalis cardinalis) established across the main Hawaiian Islands and, using this history, evaluate the role of founder effects in producing previously observed morphological divergence of these populations (Mathys & Lockwood, 2011). In the process, we also elucidate the degree to which cardinals on Hawaii have diverged from their native source population(s), and provide insight into their postestablishment population dynamics.
From written records, we know that between 1929 and 1931, 300-350 northern cardinals were purposefully transported and released onto Hawaii (Pyle & Pyle, 2009). These cardinals were shipped from the port of San Francisco (USA) and released onto Kauai, Oahu, and Hawaii Island (Pyle & Pyle, 2009). Northern cardinals are native to North America, with populations spanning the eastern half of the continent through to New Mexico and down into Mexico (Figure 1).
There are six mitochondrial clades present in North America ( Figure 1), with considerable morphological differences between them (Smith et al., 2011). The closest native population of northern cardinals to San Francisco is over 600 km to the south representing the C.c. igneous clade. There are no written records telling us whether the cardinals shipped from San Francisco came from this clade, or another one located further away but perhaps more connected to the city via train or other transportation mechanisms typical of this era. Thus, we do not know whether the cardinals throughout Hawaii were derived from one or more source clades; and, if more than one clade was involved, if a single or multiple clades founded the exotic populations on each island. The records also do not tell us how the 300-350 individual cardinals were divided across release events or how (or if) they were divided between shipments across years.
What we do know is this was the only introduction of northern cardinals to the archipelago, and they rapidly increased in population size after establishment, eventually colonizing all of the main Hawaiian Islands by the 1950s. We also know that current island cardinal populations are statistically different from each other in several morphological traits (e.g., wing and bill sizes- Mathys & Lockwood, 2011). These morphological traits are known to be heritable among birds (Badyaev & Martin, 2000a,b;Jensen et al., 2003), and Mathys and Lockwood (2011) show that observed across island differences are of such a magnitude that in situ genetic drift is not a likely causal mechanism (Mathys & Lockwood, 2011  Miles by individuals from genetically and morphologically distinct native source clades, and the morphological differences observed today recapitulate these across-clade differences ( Figure S1-Scenario 1).
Second, one or more island populations may represent an admixture of individuals sourced from different native cardinal clades ( Figure   S1-Scenario 3). Any observed differences across islands today thus evolved in response to island-specific selective forces enabled by the increases in genetic diversity that accompany admixture. Third, the cardinals on Hawaii may have been derived from one native source clade, which would suggest that current morphological differences arose after establishment from the existing genetic variation found within these founders ( Figure S1-Scenario 2).
We examine these possibilities by updating and expanding the between-island morphological trait analysis from Mathys and Lockwood (2011). We then, for the first time, compare the distribution of traits across the Hawaiian Islands to traits typical of cardinal clades in the native range. Finally, we determine the most likely native source population(s) for the exotic island populations and deduce their postestablishment population dynamics using phylogeographic and Approximate Bayesian Computation (ABC) analyses.
By combining these analyses, we assess which of the above three introduction and differentiation scenarios most likely occurred among northern cardinals in Hawaii, and we shed light on the postestablishment evolutionary dynamics of this species.

| Morphological analysis
In this analysis, we sought to establish the magnitude and direction of morphological differences in cardinals between the five main Hawaiian Islands, between the cardinals associated with each native range clade, and between the native clades and Hawaii. We used the following morphological traits: tail length, wing chord, culmen length, bill depth (at anterior margin of nares), and bill width (also at anterior margin of nares)-all measured in millimetres. These traits are commonly used metrics for evaluating evolutionary divergence between bird populations due to their known associations with life history and foraging adaptations (Lockwood, Moulton, & Anderson, 1993;Ricklefs & Travis, 1980). We visited Kauai, Oahu, Maui, and Hawaii Island in the summer of 2008, and again in the summer of 2013, to obtain morphological measurements of 74 live-caught northern cardinals. Mist nets were placed in areas that experience regular bird activity. No lures or baits were used in order to prevent bias in the sex ratio of captured individuals. Captured individuals were fitted with USGS numbered bands before release, allowing us to avoid measuring the traits of any one individual multiple times. All morphological measurements on field-captured individuals were taken in the same season, thus avoiding systematic bias in morphological traits that vary with season (e.g., wing chord- Arendt & Faaborg, 1989). Only adults were Data from live-captured individuals and museum specimens were combined for all morphological analyses. We measured only museum specimens that were captured at the same time of year as the live-caught individuals to reduce any systematic bias between the two data sources, and combined measurements for males and females to maximize sample sizes. Northern cardinals show very little differences between sexes in the traits we measured; however, to ensure that across-population comparisons were not biased by sex-specific differences, we kept sex ratios across clades and islands as close to 50:50 as possible.
Finally, it is well documented that bird specimens experience changes in some mensural characters (e.g., wing chord) after museum preparation (Bjordal, 1983;Haftorn, 1982;Winker, 1993) due to drying of the skin. In order to combine the measurements from live individuals with museum specimens, we multiplied field (live-caught) measurements of tail length and wing chord by taxon and characterspecific correction factors following Winker (1993) and Mathys and Lockwood (2011). In order to correct for individuals with missing measurements due to condition of the specimen or inability to take all measurements in the field, we approximated missing trait values using the data imputation MICE package in R (Buuren & Groothuis-Oudshoorn, 2011). This method was preferred as it has little impact on the observed population mean, uses the dataset itself to generate imputed data values, and does not reduce the variation in the dataset. We imputed trait information for <2% of the full dataset.
Recognizing that morphological traits are often intercorrelated, we collapsed the six measured traits from live-caught and museum cardinals into two principal components using Principal Component Analysis (PCA; Lockwood et al., 1993) in R statistical software with the factoextra package (Kassambara & Mundt, 2016;Team, 2014).
Prior to conducting the PCA, we log-transformed all variables and then centered and scaled the means. The first two dimensions of the PCA (PC1 and PC2) explained 75% of the observed variation in morphological traits, with PC1 capturing overall size of individuals and PC2 reflecting the ratio of the bill to body size (Table S1). We retained the PC1 and PC2 scores for each measured individual so we could compare morphological differences across populations and clades.
We initially updated and expanded the between-island morphological analysis of Mathys and Lockwood (2011) by increasing the number of individuals measured across islands and adding specimens from Maui to the comparisons. Using individual PC1 and PC2 scores, we evaluated differences in cardinal morphology between islands using multivariate analysis of variance (MANOVA) in R. If the overall MANOVA resulted in statistical significance, we followed that test with a series of pairwise MANOVA tests between islands.
Next, we compared the morphologies of Hawaiian cardinals to the four native cardinal clades for which we had sufficient data. To aid in visualizing quantitative differences in morphology across island populations and native clades, we plotted PC1 and PC2 for all measured cardinals in two-dimensional space. We visually identified individuals from each native range clade using color-coding, adding an ellipse that contained 95% of all individuals from these clades to clearly identify the range of morphologies present within each. We designated cardinals from Hawaii with a unique color code as well as designated individuals according to their island of residence using island-specific symbols. This graph allows one to visualize the morphological "map" of native range cardinals, where each clade occupies a relatively distinct position in the two-dimensional space, and then visually evaluate where the Hawaiian cardinals "fit" onto this map.
Using these data, we evaluated the following scenarios: (1) the Hawaiian cardinals fall entirely into the trait space of only one native range clade, indicating all Hawaiian cardinals were derived from this single native source and any divergence they show across islands is typical of the range of morphologies seen in that clade; (2) Hawaiian cardinals span two (or more) native range clade spaces and that cardinals from one island clearly fall within one clade and cardinals from another island fall in the other clade, indicating that interclade morphological differences are being recapitulated across islands (founder effect); or (3) the Hawaiian cardinals do not neatly fit into any single native clade's morphological space, indicating potential admixture at the time of introduction, postestablishment divergence, or both. We quantitatively evaluated differences in PC1 and PC2 between clades and Hawaii with MANOVA followed by pairwise MANOVA.

| Sequence data generation
In order to determine the native source(s) of cardinal populations across Hawaii, we combined the C. cardinalis native range genetic data from Smith et al. (2011) and genetic information from the live-caught individuals to create a merged northern cardinal dataset. Smith et al. (2011) used the sodium dehydrogenase subunit-2 (ND2) mitochondrial locus to establish discrete genetic boundaries for six native range clades. In order to compare Hawaii cardinals with this dataset, we used the same locus. We collected feathers from 46 of the measured individuals caught in the Hawaiian Islands in 2013.
Feathers were placed in small envelopes, and upon return to the lab, the calamus from multiple feathers was clipped to obtain a biological sample for each individual. These samples were placed in 1.5ml Eppendorf tubes in order to extract genomic material from cells found on the feather calamus. We extracted DNA using a DNeasy blood and tissue kit under standard protocols (Qiagen reference), with Proteinase K incubation taking place overnight (minimum of 8 hr) to ensure complete digestion.
We amplified 1,042 base pairs of ND2 via polymerase chain reaction Piscataway, NJ, USA). Sequences were obtained using both primers to create a consensus of the full 1,042 bp ND2 sequence after chromatograms were cleaned and aligned in Sequencer 5.1 (GeneCodes, Ann Arbor, MI, USA). All sequences were evaluated for insertions and deletions, as well as translated to amino acids to check for stop codons and the presence of nuclear copies (Sorenson & Fleischer, 1996).

| Phylogeographic analysis
We executed a phylogeographic analysis using the merged northern cardinal dataset to determine which native source clades were associated with each exotic island population. We ran the dataset through the program PartitionFinder 1.1.1 (Lanfear, Calcott, Ho, & Guindon, 2012) in Python v2.7 under two model schemes: unpartitioned whole gene ND2 sequences and partitioned by codon.
We implemented a MrBayes model filter to select only the twentyfour DNA evolutionary models that were compatible with the MrBayes program. PartitionFinder generated model schemes for both partitioned and unpartitioned data and ranked them by Akaike Information Criterion (AIC). We then constructed a phylogenetic tree with MrBayes v3.2.2 (Ronquist & Huelsenbeck, 2003) under the selected best scheme for both unpartitioned and partitioned data. The program was allowed to run for 10 million generations, while being sampled every 1,000, with a relative burn-in of 0.25. We visually inspected MCMC chains using the program Tracer v1.6 (Rambaut, Suchard, Xie, & Drummond, 2014) to confirm adequate burn-in and convergence of chains and used FigTree v1.4.2 (http://tree.bio.ed.ac. uk/software/figtree) for final tree assembly and inspection.

| Approximate Bayesian computation analysis
We used Approximate Bayesian Computation (ABC) to test a suite of possible introduction and range expansion scenarios. Briefly, ABC is a Bayesian analysis that allows for direct comparison of multiple introduction hypotheses (known as scenarios) and provides relative probabilities for each, given the data provided. This comparison is accomplished by performing inference computations from simulated pseudo-observed datasets (PODs) that take into consideration the putative introduction histories modeled, moving backward through time from the observed data. The PODs most similar to the observed dataset are then selected (with replacement) via a Euclidean distance measure (Cornuet et al., 2014;Lombaert et al., 2010;Valentin, Nielsen, Wiman, Lee, & Fonseca, 2017). The selected PODs have relative posterior probabilities calculated for their respective scenarios via a logistic regression estimate, allowing the user to select a significantly different scenario as being most likely to have occurred (Cornuet et al., 2014;Valentin et al., 2017).
We framed testable scenarios around three main questions: (1) Can we identify from which of the source clade(s) the cardinals brought to Hawaii (the founding cardinals) were sourced? (2) Can we assess if the 300+ cardinals that reached Hawaii were effectively divided into three evenly distributed groups of founders and released simultaneously across all three islands, or were approximately 100 founders introduced to a single island during each introduction event over 3 years? and (3) Can we identify which island population(s) provided the founders of the Maui population? For each question, we modeled two or more scenarios and then compared these against each other in order to quantify their relative probabilities.
We used the program DIYABC to conduct these analyses (Cornuet et al., 2008(Cornuet et al., , 2014 and used the following summary statistics to conduct our analyses: one-sample statistics-number of haplotypes and number of segregating sites, two-sample statistics-number of haplotypes. To address the first question (source clade), we evaluated four variations of three scenarios. The first scenario supposed that the source of Hawaii cardinals was the western region of the source clade (see Results for clade analysis below; Figure 1). The second scenario supposed the source individuals were derived from the eastern region of the source clade (Figure 1). The third scenario supposed that the Hawaii population was a mix of both regions.
For these three scenarios, the first variation evaluated which region was the likely source of the Hawaii introduction without enforcing a change in effective population size (i.e., no genetic bottleneck). The second variation reduced the effective population size after initial introduction into the Hawaiian archipelago (genetic bottleneckconditioned to be less than both native sources) and then allowed the population to change (no condition set to Hawaiian populations).
The last two variations (three and four) considered the possibility that each source region contained an unsampled population that was the source of Hawaii founders, and contains genetic haplotypes not present in our dataset. Variations three and four were identical to the above second and first variations, respectively, except an unsampled population for each region was used rather than the region data itself. The variation with the highest confidence in scenario choice (i.e., contained the least amount of error) and contained a statistically significant scenario was considered the most probable, given our data.
To address the second question (pattern of release events), we evaluated two variations of two scenarios. The first scenario supposed the 300+ cardinals transported from the mainland were equally divided among the three islands, but equal subsets were released in 1929, 1930, and 1931 resulting in smaller founding population sizes. The second scenario supposed that of the 300+ founding cardinals, roughly 100 were acquired and introduced to one island per year. We again evaluated whether there was evidence of a population bottleneck with our scenario variations. For the first variation, there was no change in effective population size enforced after founders were introduced to Hawaii (i.e., no enforced bottleneckno restrictions placed on Hawaii parameters). For the second variation, we did enforce an initial reduction in effective population size (i.e., bottleneck-restricted Hawaii parameters to be less than native range and fit scenario) and then allowed the population to increase.
To address the third question (source of Maui cardinals), we evaluated three different scenarios: (1) colonizers to Maui came from Hawaii Island; (2) colonizers came from Oahu; and (3) colonizers were derived from both islands.
In all ABC scenarios, we set parameter priors to fit a uniform distribution (under default bounds) and placed conditions on parameter priors only to fit the intention of each scenario as defined above.
We chose the HKY mutation model, based on the results from PartitionFinder during the phylogenetic analysis (see Results below, Table S2), and set it identically for all scenarios evaluated (Table 1).
We ran all experiments for three million computations prior to conducting any analyses.

| Morphological analyses
Reinforcing the findings of Mathys and Lockwood (2011), we found that northern cardinal populations showed substantial morphological divergence across the main Hawaiian Islands (Table 2, Figure 2).
In particular, cardinals from Hawaii Island differ from those on all other islands except Maui (Table 2). Cardinals resident on Hawaii Island and Maui tend to be larger than their counterparts on Oahu and Kauai, especially in tail length (Figure 2). We also find residents of Maui have significantly larger wings than all other Hawaiian island populations ( Figure 2).
Our evaluation of morphological differences among the four evaluated phylogeographic clades confirms the existence of substantial morphological variation between northern cardinal clades across their native range (Table 2, Figure 3). In particular, we found that the native populations differ substantially in body size with C.c. igneous being the largest of the set, C.c. cardinalis moderately large sized, and the two lower-latitude clades in Mexico the smallest (Figure 3). Bill dimensions also vary across clades, with cardinals exhibiting somewhat shorter and pointier bills (relative to body size) in the southern Mexican clades as compared to the two clades that cover sections of the United States  for which we were able to obtain measurements. In addition, there is no clear pattern whereby the morphology of cardinals resident on an island corresponds to the morphology of cardinals from any one clade (Figure 3). Thus, there is no indication from this analysis that the pattern of morphological divergence observed on the islands matches any observed pattern of morphological differentiation among clades across the native range.

| Sequence generation and phylogeographical analysis
After amplifying and sequencing the ND2 mtDNA locus for the 41 northern cardinal samples obtained from Hawaii, we found a total of 19 haplotypes (Table 3). We observed six, ten, and seven haplotypes for the cardinals present on Oahu, Hawaii Island, and Kauai, respectively, while Maui had just three haplotypes. These sequence data can be found in Genbank under accession numbers MH010209-MH010303.
After combining our sequence data with that of Smith et al.   (Figure 4), which also did not show geographical assortment across the wide sampled range ( Figure S2). F I G U R E 3 Two-dimensional representation of northern cardinal morphological trait space using PC1 and PC2. PC1 reflects overall body size, whereas PC2 measures how bill depth and width change relative to a change in body size. We only include four of the native range clades in this analysis due to low sample size in two clades. Each oval encapsulates 95% of the variation in morphology between the individuals we measured, representing a clade-specific trait space. Large symbols within each oval depict the mean PC scores for each clade. We depict all individual cardinals captured and measured in Hawaii in light blue (with different shapes for each island) to visually show the distribution of their morphology (light blue oval) relative to native clades PC2 PC1

| Approximate Bayesian computation analysis
The first question was aimed at identifying the region within the native source clade (C.c. cardinalis) from which the Hawaiian cardinals were derived. Without an enforced bottleneck, we did not find a significant difference in the relative probabilities among any of the three tested scenarios (Table 4). When we enforced a population bottleneck (variation two), however, we found the scenario where cardinals on Hawaii were derived from the eastern region of the C.c. cardinalis clade to be significantly more likely (Table 4). However, both of these variations had confidence scores below the remaining two, which included unsampled populations from each region. Of the remaining two variations, the third had the highest confidence score (0.643, Table 4

| D ISCUSS I ON
A species' invasion history can profoundly influence the degree of divergence observed, including via founder effects whereby phenotypic and genetic spatial structure in the species' native range is TA B L E 3 Summary of number of samples (n.) used in the genetic analyses conducted, with localities sorted by mtDNA clade for the native range (with the west and east regions for C.c. cardinalis identified) and the Hawaiian archipelago. Each clade, and Hawaii, is further subdivided by locality, while providing the number of haplotypes per location (n. Haps) and the haplotypes observed. Any haplotypes followed by a number in parentheses indicates multiple specimens observed with said haplotype.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N
REV performed all laboratory work, collected morphological data, performed all morphological and genetic analyses, and wrote the manuscript. JLL collected morphological data, oversaw all morphological analyses, and wrote the manuscript. BAM collected morphological data, collected biological material for laboratory work, and wrote the manuscript. DMF oversaw all laboratory work, oversaw all genetic analyses, and wrote the manuscript.