Origin identification of migratory pests (European Starling) using geochemical fingerprinting

Study Area

British Columbia (BC) is the western-most province of Canada and has high geological diversity in rock types (Natural Resources Canada, 2012) (Fig. 1). The Okanagan-Similkameen region in south-central British Columbia is about 200 km long and 20 km wide (Fig. 1). It lies in between the Columbia and Cascade Mountain ranges in the southern interior of BC. The valley bedrock is comprised of volcanic, sedimentary, metamorphic and intrusive rocks (Okulitch, 2013). Southern Okanagan-Similkameen region is heterogeneous, made up of volcanic, metamorphic, plutonic, and sedimentary rocks, and overlain by glacial sediments (Okulitch, 2013). Central Okanagan is made up lava, andesite with quartz-filled amygdales. The northern Okanagan is made up of Paleozoic phyllite, argillite, schist, amphibolitic, quartzite, carbonaceous, micaceous chlorite and other minor constituents (Table S1) (Okulitch, 2013).

Starling habitat is generally the lower elevation portion of a long valley, approximately between 300 and 700 m elevation. Above 700 m, the valley rises quickly to a plateau of >1,200 m (Fig. S1), greatly constraining distribution and likely movement and migration to a relatively small area (Fig. S1).

Sampling

Sampling was performed in summer (May–July) and in fall/winter of 2015. In summer, only juveniles, having little prior chance to disperse beyond their natal habitat, were collected to establish fingerprints of source populations. Thus, their bone tissue will reflect the elemental fingerprint of the sites where they were collected. To make our trace element fingerprint library more representative, juveniles were collected from almost all the Okanagan-Similkameen cities/towns, and lands used for different agricultural activities, such as vineyards, orchards, dairy farms, feedlots, grasslands, and landfill sites.

In fall, problem birds (both adult and juveniles (<1 year)) were collected in vineyard/orchards and dairy farms/feedlots. The age of the birds was determined from their plumage (Szép et al., 2003; Szép et al., 2009). Additionally, aging by plumage was validated by the chronology of skull development (skull ossification–skulling) in juveniles (McKinney, 2004). Five randomly selected samples, which were considered as juvenile through plumage, were dissected to observe developmental stage by skulling (Mueller & Weise, 1996). These starlings are considered pests and hereafter referred to as “problem birds”.

Sampling was performed in collaboration with the Starling Control Program (SCP), British Columbia Grapesgrowers’ Association (BCGA). The work was conducted in accordance with the ethics training requirements of the Canadian Council on Animal Care (CCAC)/ National Institute Animal User Training (NIAUT) program certificate number 7252-1. A team of professional trappers contracted by the BCGA, has been trapping and euthanizing starlings following the guidelines of the Canadian Council on Animal Care every year. We obtained frozen sample from SCP trappers and each individual bird is considered as a single sample. Altogether, 105 juveniles were collected from ten different locations situated at distances from 12 to 190 km apart throughout the Okanagan-Similkameen region in 2015. The number of juvenile samples collected in summer varied depending on availability at each location, 20 from Kelowna, 10 from Hullcar, 10 from Salmon Arm, 9 from Armstrong, 9 from Mara, 18 from Vernon, 10 from Oliver, 8 from Osoyoos, 6 from Penticton and 5 from Keremeos. In fall, 118 problem birds (including young adults <1 year, n = 20; and adults: 1+ years, n = 98) were collected from vineyards/orchards and dairy farms and feedlots having major starling problem. Some of the problem bird sampling locations are the same where juveniles were collected like Vernon (n = 4), Pentiction (n = 27), Keremeos (n = 5), Kelowna (n = 14), and Oliver (n = 15); however, Okanagan falls (n = 4), Summerland (n = 30), Coldstream (n = 10), and Lumby (n = 9) were added in fall. These locations spread across Okanagan-Similkameen region from south to north (Fig. 1).

Sample preparation

Samples were prepared by methods described by Norris et al. (2007) with some modification. Briefly, the bone was cleaned by removing outer skin and marrow. Bones were washed repeatedly with ultra-pure water (18 M Ω cm⁻¹), dried at room temperature (20 °C) (Szép et al., 2009) and weighed. The clean fragments of the bone sample, weighing between 30 mg to 80 mg each, were put separately into Teflon tubes. Trace element grade concentrated nitric acid (2 ml) was added to each tube and placed in a heating block 75 °C (±5 °C) until the bone was completely dissolved. Samples were then cooled to room temperature, and then one milliliter of trace metal grade hydrogen peroxide was added to further digest the dissolved organic compounds (Norris et al., 2007). Samples were then evaporated on a hot plate at 85 °C (±5 °C) (Szép et al., 2003; Neuhauser, 2013). After cooling, 2 ml of 1% HNO₃ with 1 ppb indium internal standard was added to each vial to dissolve the residue over a period of 2 h at room temperature.

Samples were transferred quantitatively with three washes of the HNO₃ indium solution into acid-washed high-density polyethylene tubes to a final volume of 10 mL (Norris et al., 2007; Neuhauser, 2013). Sample blanks were prepared in the same way for every set of digestions to monitor contamination (Donovan et al., 2006; Norris et al., 2007; Szép et al., 2009). Sub-samples were further diluted 100-fold with 1 ppb indium-1% HNO₃ for ICP-MS analysis and 5 and 10 fold with 1ppm yttrium-1% HNO₃ for ICP-OES analysis. Samples were diluted to fall within the detection limits of the instruments.

Trace element analysis

We measured 21 elements: aluminium (Al), silver (Ag), barium (Ba), calcium (Ca), cadmium (Cd), cobalt (Co), chromium (Cr), copper (Cu), molybdenum (Mo), manganese (Mn), lead (Pb), sulphur (S), scandium (Sc), selenium (Se), tin (Sn), strontium (Sr), vanadium (V), zinc (Zn), magnesium (Mg), sodium (Na), and potassium (K) in bone tissue. The elements were selected on the basis that they are present at levels high enough for reliable quantification, and could also vary in space.

Seventeen of the 21 elements were analysed by Thermo-Fisher Element XR sector field Inductively Coupled Plasma-Mass Spectrometer (ICP-MS) (Poesel et al., 2008; Neuhauser, 2013). A subset of four highly-abundant elements, Ca, Mn, Na, and K, were analysed via a Thermo-Electron Corporation, iCAP 6000 Series XR Inductively Coupled Plasma-Optical Emission Spectrometer (ICP-OES) (Szép et al., 2003; Szép et al., 2009; Neuhauser, 2013). For ICP-MS, elements were analysed in either low resolution or medium resolution mode to resolve polyatomic interferences. Both instruments were calibrated by external multi-element standards. Four different concentrations of multi-element standards were run to produce calibration curves (Poesel et al., 2008). Indium (Norris et al., 2007) and yttrium (Neuhauser, 2013) were used as an internal standard for ICP-MS and ICP-OES respectively to correct for instrument drift, sample density and to improve accuracy and repeatability. The blank and reference solutions were run at 20-sample intervals throughout the analysis. Three replicate analyses were done for each sample. Relative Standard Deviation (RSD) of the analysis was between 1–2%.

Statistical analysis of data

The elemental concentration was normalized to the concentration of calcium rather than mass to minimize variability in the ratio of mineral:collagen content. The collagen structure of bone bonds relatively poorly with metals, whereas the mineral apatite substitutes cations and anions relatively easily into its structure. As birds develop, the relative content of apatite increases (Rogers & Zioupos, 1999; Currey, 2004; Currey, Brear & Zioupos, 2004; Pasteris, Wopenka & Valsami-Jones, 2008). Normalizing with calcium, therefore, removes this variability due to collagen content. After calcium normalization, the data were then standardized via z-scoring to prevent individual elements from having a disproportionate influence on the groupings derived from multivariate analyses (Fowler, Cohen & Jarvis, 1998; Norusis, 2016).

A multivariate analysis of variance (MANOVA) was used to analyse trace element profiles of juvenile bone tissue collected in different locations (Norusis, 2016). A multiple discriminant analysis (Hair et al., 2010) of juvenile samples was performed to discard elements (variables), which were little related to group distinction, and also to develop the predictive model of group membership based on trace elements (chemical profile of bone). A discriminant function (equation) was developed based on the linear combination of the predictor variables (trace elements) that provides the best discrimination (correctly separating individual birds) among the predetermined groups (Lachenbruch & Goldstein, 1979; Greenough, Longerich & Jackson, 1997; Szép et al., 2009). The predetermined groups were the sites where the juvenile starlings were caught. The stepwise method and Mahalanobis distance (Hair et al., 2010) were used for the analysis because of the large number of variables (20 different element concentrations).

The stepwise method looks at each element (variable), one at a time, and determines which element is the best predictor of group membership, and ultimately generates the best set of variables to predict group membership. Cross-validation assessed the success of the proper bird groupings via the discriminant function and rules. In order to classify objectively the problem birds with respect to origin, the value obtained from each function for each problem bird was normalized (ratio of sample/ratio of centroid) by the centroid value of each site. Thus, if an individual problem bird were to lie at the centroid for a particular site, its value with respect to that site would be one. Each bird’s probable origin was identified by how closely that bird’s normalized value approximated the centroid values of each site. For some of the problem birds, the discriminant function value was too high to assign a specific juvenile source population. For such birds, we have categorized them as immigrant/unassigned birds.

Cluster analysis (Hartigan, 1975) was performed to evaluate the spatial separation of the elemental fingerprint of starlings and the identification of likely immigrants. In this study, the ward.D method (Ward, 1963; Murtagh & Legendre, 2014) with Euclidean distance measures was used for analysis, which provided better separation of individual birds with similar signatures and geographic location. Additionally, an average method with correlation distance was used to compute a test of significance of the clusters (by calculating the p-value of the cluster). The R package “pvclust” (Suzuki & Shimodaira, 2014) was used; it uses bootstrap resampling techniques to compute the p-value for each cluster. This method generates thousands of bootstrap replications by randomly sampling elements of the data. For each of the clusters, Approximately Unbiased (AU) and Bootstrap Probability (BP) values were calculated. The AU probability values (p-values) are computed by multiscale bootstrap resampling where AU ≥ 95% are considered to be strongly supported by the data, while BP corresponds to the frequency that the cluster is identified in bootstrap copies (Suzuki & Shimodaira, 2014).

Statistical analysis was conducted in SPSS ver. 24.0 and R statistical software version 3.3.1 (R Development Core Team, 2014) with the cluster (Maechler et al., 2017) and pvclust (Suzuki & Shimodaira, 2014) packages. Map projections of data were done in ArcGIS 10.4.

Spatial separation of source population of starlings in the Okanagan-Similkameen region

The source population trace element fingerprint library for the Okanagan-Similkameen region was developed from juvenile starlings collected from 10 different locations. Multivariate comparison of 20 elements in the juvenile (n = 105) population using MANOVA shows there is a significant difference in the trace elemental composition of juvenile bone among locations (Wilks’ λ <0.0001, F = 6.655, df = 180, 645.230, power = 1, p < 0.0001).

A stepwise discriminant analysis of 20 trace elements based on the juveniles (n = 105) collected at 10 different locations, provided four canonical discriminant functions with significantly high eigenvalues (>1), and canonical correlations (rc > 0.70), that separated juvenile samples from respective sites (Wilks’ λ = 0.001, 0.11, 0.73, χ² = 633, 428, 247, p < 0.0001) (Fig. 2). Ten elements (K, Mg, Na, Ag, Cd, Ba, Sc, Cr, Cu, and Se) out of 20, were used as the best predictor variables in discriminating group membership. The probabilities of correctly classifying juveniles to their respective locations using stepwise discriminant analysis are shown in Table 1. The table shows the predicted classification (group membership) using cross-validation from different locations compared to the apparent or actual observation. The overall probability of correctly classifying juveniles was 79% (P < 0.0001, Press’s Q; Table 1). The juveniles were grouped 67%–100% correctly, with Osoyoos as an exception. Only 25% of the Osoyoos juveniles were correctly grouped by location. Mostly the misclassification pattern that occurred was among nearby locations. For example, a high rate of misclassification, 30% and 37% (Table 1), occurred between Oliver and Osoyoos (∼17 km apart), 33% between Penticton and Keremeos (∼45 km apart), and 22% and 33% between Mara and Armstrong (Table 1). Hence, almost all the misclassification occurs between sites that are very close to each other, with one exception, between Osoyoos and Vernon in 2015 (37%, ∼165 km apart, Table 1).

Cluster analysis of juveniles shows broadly similar results. The clusters (ward.D method) grouped juveniles from different sites into different clusters, even separating individuals derived from locations ≤12 km (Fig. 3) apart. Since there is no particular rule for choosing the height in the ward.D method dendrogram, the rule of thumb height of 10 was applied. Between 33–100% of juveniles from the same locations were grouped together at the height of 10 (Fig. 3). Approximately unbiased (AU) values obtained via the average method were used to test the significance of the clusters. Ten major significant clusters with red boxes (AU p-value > 0.95) separated juveniles from different geographic locations (Fig. S2). Thus, juvenile grouping with respect to geographic location in cluster analysis supports the predicted membership results derived from stepwise discriminant analysis. Thus, our results are robust to different multivariate statistical methods.

Identification of the source population of problem birds in the Okanagan-Similkameen region

The source population of problem birds was identified by comparing their trace element fingerprint with the juvenile fingerprint library. Cluster analysis of all birds including both juvenile (n = 105), and problem birds (<1 year adults (n = 20), and 1 + year adults (n = 98)) depicts two distinct clusters, separating most problem birds into one cluster and another cluster comprising a mix of the remaining problem birds and juvenile birds (Fig. 4). Around half of the problem birds (n = 52, 44%) clustered separately, indicating chemical composition distinct from the juvenile fingerprint library (Fig. 4). Problem birds having trace element fingerprints distinct from the Okanagan-Similkameen populations (doesn’t match our library) were considered as immigrant/unassigned birds, likely originating from unsampled sites within the Valley or from outside of the Valley.

To identify the origin of the remaining (n = 66, 56%) local problem birds in the Okanagan-Similkameen region, the same four canonical discriminant functions were used. The first function explained 44.2% of the total variance among locations. The second function explained 33.0%, and the third and fourth explained 11.7% and 6.1%, respectively of the total variance among locations. The problem birds collected in fall were divided into two groups. First, problem birds collected in vineyards and orchards and second, those collected at dairy farms and feed lots. The majority (55%) of problem birds caught in vineyards and orchards were consistent with trace element signatures of juvenile birds from the North Okanagan. Specifically, the highest contributing site was Vernon (41%) (Fig. 5). The remaining 45% of problem birds were a mixture of local and immigrant/unassigned birds (29% southern local, 2% central local), 14% of problem birds in vineyards and orchards were immigrants/unassigned (Fig. 5). In contrast to vineyards and orchards, the problem birds caught on dairy farms and feedlots were largely immigrants/unassigned (89%) and only 11% were local from northern site (Vernon) (Fig. 5).

Spatial separation of starlings

Resolving populations in space revealed distinct trace elemental signatures within the 200 km length of the Okanagan-Similkameen region. The spatial separation of starling in this study using bone tissue trace elemental fingerprinting is comparable to spatial separation of other birds at scales ranging between kilometers to continental distances using feather tissue trace elemental fingerprinting (Poesel et al., 2008; Szép et al., 2009). For example, white-crowned sparrow (Zonotrichia leucophrys F.) samples collected from four dialect populations were distinguished within 400 km (Poesel et al., 2008). Similarly, sand martin (Riparia riparia L.) and barn swallow (Hirundo rustica L.) populations were separated within and between continents (i.e., across Europe and between European and African moulting sites) (Szép et al., 2003; Szép et al., 2009). In our study, the shortest distance separating sites with distinct chemical signatures was 12 km (between Vernon and Armstrong).

Site-specificity of the elemental compositions of bones from juvenile starling is similar to analyses of natal origin from nestlings and juvenile sand martins (Riparia riparia L.) (Szép et al., 2003) and for western sandpiper using feather tissue (Calidris mauri Cabanis; Norris et al., 2007).

Moreover, a high degree of similarity in feeding and foraging habitat of juveniles could be factors maintaining the similarity of juvenile bone trace element composition in local populations. Since most juveniles collected within sites were chemically similar, it appears that these juveniles had not dispersed. Thus, the juvenile fingerprint library adequately represents the fingerprints of source populations at particular sites in the valley.

The chemical signature of a few juvenile birds caught within a site matched with fingerprints from other sites according to the predicted group membership derived from discriminant analysis. Such misclassification of juveniles among nearby sites such as Oliver and Osoyoos may have arisen due to overlapping foraging ranges of the juveniles or adults feeding chicks, geological similarity between nearby sites, similar agricultural practices and food choice. These mechanisms are not mutually exclusive and all may have been operating. Finally, misclassification noted in a few birds between the sites that are 100s of kilometers apart geographically might be due to differences in food choice and early short distance dispersal or geochemical similarity between remote sites.

The overall degree of classification of juvenile fingerprints was 79% among sites. This is similar to the degree of differentiation reported for shearwaters (Calonectris Cory & Oustalet) 75-89.9% (Gómez-Díaz & González-Solís, 2007) and lower than that in white-crowned sparrows 100% across 400 km (Poesel et al., 2008). The dissimilarity of chemical composition among sites within the Okanagan-Similkameen region may be due to geologic factors (Ethier et al., 2013), and/or hydrologic factors including irrigation (Poesel et al., 2008). The underlying geology of the Okanagan-Similkameen watershed is very diverse within small spatial scales (Okulitch, 2013).

Similarly, soil chemistry is non-uniform over the landscape (Goitom Asfaha et al., 2011). Water, mineralogy and soil microbiology affect soil chemistry. Solutes present in soil water can be derived from weathering of rocks and minerals locally or could have come from water (precipitation and irrigation). The trace elements move from minerals to soil water (Galgano et al., 2008) where they are taken up by plants. In turn, plants form the base of the food web, eventually accumulating in animal tissues, forming a unique marker. The trace element availability in soil can be further modified by factors like soil pH, water content and porosity (Kelly, Heaton & Hoogewerff, 2005; Kim & Thornton, 1993).

Trace element composition in bone tissues can yield a very specific biogeochemical marker. Out of 21 elements analysed, ten elements (K, Mg, Na, Ag, Cd, Ba, Sc, Cr, Cu, and Se) were used in a discriminant analysis that classified juveniles source populations. All are cations (Se is also an oxyanion) and could potentially substitute into apatite. For example, Pb substitutes for Ca in hydroxyapatite in bone with a turnover rate ranging from 1–8% per year depending on the type of bone tissue (Martin et al., 2007). The shortlisted elements from the discriminant analysis are therefore a mixture of those elements, which could easily be present in bone, and could also vary in space.

Assigning problem starlings to natal populations

Problem birds were assigned to potential natal populations through discriminant analyses. However, some of the problem birds could not be assigned to its source population and were designated as immigrant/unassigned problem birds. Failure to assign origins for a portion of the problem birds could be caused by continued addition of elements to bone over large unknown space over time. Second, adult starlings >1 year old, natal fingerprints may differ inter annually, especially if foraging by such a generalist might change for different climatic years. Third, problem birds could come from outside the region.

The problem bird contribution to vineyards and orchards in the southern part of the Okanagan-Similkameen region was mostly from the northern end of the valley. It is likely that the birds reared in the northern part of the valley move south as fruit ripens in the fall. The remaining problem birds in vineyards and orchards are a mixture of local (southern) and immigrant/unassigned birds. This is consistent with documented migration patterns in starling. Generally, fall migration occurs from September to December (Kessel, 1953; Dolbeer, 1982).

The problem birds on dairy farms and feedlots were a mixture of immigrants/unassigned and local birds. The high immigrant/unassigned population in the north Okanagan is likely because these sites are close to unsampled source populations further to the north. During winter, dairy farms and feedlots provide, adequate shelter and warmth and a diverse supply of food and water (Palmer, 1976). Several past studies have noted that the flocks of starlings increase in number and concentrate on dairy farms and feedlots in winter (Palmer, 1976; Glahn, 1981). Moreover, one of the challenges to this type of studies is that the trace element profiling is more difficult for wide-ranging, continuously distributed species as the number of potential source areas is vast. If the organism migrates long distances, then trace elements alone might not be sufficient to track their origin. Also trace elements alone would not be sufficient for tracking the origin of migratory organisms that undertake large-scale movements and forage in diverse geochemical environments. However, the combined use of trace elements and isotopes could help track continental-scale movement patterns. Besides identifying bird origins, the trace elemental signature could be used in several other ecological studies of birds, including matching breeding populations with overwintering populations, tracking fine-scale dispersal (Poesel et al., 2008), identifying migratory connectivity (Szép et al., 2009) and identifying moulting areas (Szép et al., 2003).