Using trace elements to identify the geographic origin of migratory bats

Soil data collection

To generate a trace element basemap, data on trace element composition of soil were collected from the US Geological Survey (USGS) and Natural Resources Canada (NRC) online databases available at https://mrdata.usgs.gov/soilgeochemistry/#/summary and https://www.nrcan.gc.ca/home, respectively. We used data that had originally been collected as part of a continent-wide low-density geochemical cooperative survey between USGS, NRC, and Mexican Geological Survey (Smith et al., 2005; Smith et al., 2014) since previous studies have found little difference between sampling at low-density and high-density on soil trace element interpolation (Birke, Rauch & Stummeyer, 2015). However, we were only able to access data for US and Canada. Details of the sampling techniques are found in Smith et al. (2014). In brief, 4857 low density sites (1 site per 1,600 km²) were selected and samples collected from soil layers of 0–5 cm, Soil A (topsoil), and Soil B+C composite (deeper soil up to 1m). We used trace element concentration from Soil B+C composite samples to minimize the effects of potential surface contaminants (Sastre et al., 2001). Elements selected for analysis were the same as those chosen for analysis in fur samples (see below for inclusion criteria) and included Aluminum (Al), Nickel (Ni), Copper (Cu), Rubidium (Rb), Yttrium (Y), Molybdenum (Mo), Tin (Sn), Barium (Ba), Cesium (Cs), Cerium (Ce), Mercury (Hg), Magnesium (Mg), Manganese (Mn), and Iron (Fe).

Fur analysis

To validate the use of trace elements for sourcing individual bats it is necessary to link a resident bat’s fur elemental profile to that of a specific geographic region. To accomplish this, we used fur samples from eastern red bats collected from on axillary region of 126 museum specimens from the Smithsonian National Museum of Natural History that had previously been used for a similar purpose in an isotope validation study (Pylant, Nelson & Keller, 2014; Fig. 1). Permission to use these samples in this study was granted by the Smithsonian (J Wieringa, pers. comm., 2019). Following guidelines from isotopic studies (Cryan, Stricker & Wunder, 2014), we assume these samples are representative of the element profiles of local bats because these samples were collected during the summer, presumably after molting and prior to fall migration (Fraser, Longstaffe & Fenton, 2013). We assume that the fur collected was grown at or near the sampling location. To sample each individual, fur sample was clipped from the torso of each bat and placed inside 2 ml centrifuge tubes for storage until later use.

To prepare samples for trace element analyses, individual fur samples were processed following methods in Flache et al. (2015a) and Hickey et al. (2001). Specifically, each sample was cleaned with Triton-X (1:400) and acetone, rinsed with reagent grade water and dried at 60 °C in a sterile Low Density-Polyethylene (LDPE) tube. Gold chloride (AuCl₃) was then added at 1 ppm final concentration as a mercury stabilizer (United States Enivronmental Protection Agency, 2003). We then added 1 ml of HNO₃ and placed tubes on a heating block at 70 °C for 1 h to ensure complete fur digestion. Samples were then allowed to cool to room temperature, after which 10 ppb final concentration Indium was added as an internal standard, and samples were then diluted with reagent grade water to 10 ml.

To more accurately determine trace element concentrations, we used a narrow range calibration curve (calibrations limited to only expected range) as determined from concentrations of previous analyses (Table S1). To create our calibration curve, 5 bat fur samples were analyzed using ThermoFinnigan Element 2 High Resolution ICP-MS at the Trace Element Research Lab (TERL) at Ohio State using the Semi-Quantitative mode (Chen et al., 2008; Table S1). This mode allows for rapid quantification of all elements through estimating only a few elements accurately and inferring the rest on previously determined ratios. While there is increased error, this method of analysis allows for rapid quantification to assess which elements to use in further analysis. When compared to blanks, these 5 samples identify elements that were reliably detected above background and gave concentrations that were used to create narrow range calibration standards. This more rapid analysis allowed us to better determine which concentrations to expect for more in-depth determination (below) which created more accurate calibration curves, and allows for better determination of concentrations.

We used the normal quantification mode on the ThermoFinnigan Element 2 ICP-MS to make additional elemental concentration measurements. Samples were fed into the ICP-MS with a delay before reading to reduce the chance of carry over effects, followed by a short wash between each sample. Samples were then ionized in plasma and the concentrations determined when compared to a calibration curve. After every 10-15 samples, one of the calibration standards used to create the calibration curve and a blank were re-analyzed to ensure accuracy of our results. Additionally, ten samples were re-run to determine analytical variation. We used software included with the ICP-MS to correct each sample to the values for the internal standard and further corrected sample concentrations based on the mass of the fur used.

Statistical analyses

Key steps in our statistical analyses were as follows: (1) Summarize trace element variation among individual bat fur samples using spatially explicate measures; (2) Link fur and soil elemental concentrations across the range of this species; (3) Generate a basemap across the species range; (4) Assign individuals to origin, and assess the accuracy and precision of assignment.

Summarizing individual variation

Sourcing individuals via their trace element profiles requires summarizing the spatial correlation of trace elements among individuals sampled (Donovan et al., 2006). Therefore, we used sPCA scores to summarize variation for bat fur concentrations using the function ‘gwpca’ in the ‘GWmodel’ package in R (Gollini et al., 2013). sPCA introduces a spatially based weighting matrix into the PCA matrix. The advantage of sPCA over traditional PCA is that sPCA maximizes spatial autocorrelation, and therefore shows strong spatial relationships in the first few principle axes. When using ‘gwpca’ we generally applied the default values in the package, with the exception that bandwidth was determined using the ‘bw.gwpca’ function. This function plays an important role in the weighting function of the sPCA. The loadings were then used to determine the principle component scores based on element profiles for individual bats. These values could then be used for geographic assignment; both for creating a basemap, and validating their use.

Linking fur and soil trace element values

We first tested the assumption that element values from fur reflect values in the soil from the same geographic location. To do this we generated linear models for individual and multivariate measures of element abundance which estimated the relationship between soil element concentration (extracted from maps in Smith et al., 2014) and fur element concentration. This was accomplished using the function ‘lm’ in R (R Core Team, 2019). The usefulness in assignment was assessed by the amount of variation explained by the linear model (R²) and the significance of the relationship. However, we found weak relationships for individual elements, which precludes the use of single elements for assignment. Weak univariate relationships does not, necessarily, preclude significant multivariate relationships from being present. In which case, principle components generated above may capture patterns of variation across multiple elements that are useful for assignment. To assess the ability of use of a composite metric (principle components), we compared predicted and actual values of PCs, and compared soil to hypothetical fur to PCs (described below) to actual fur PC values.

Geographic assignment

We used methods developed for sourcing species using hydrogen isotope values of fur (e.g., Pylant et al., 2016) to conduct geographic assignment analyses. We created a basemap of multivariate estimates of fur values in ArcMap 10.6, after removing outliers using squared robust Mahalanobis distances and evaluating plotted values (Ethier, Kyle & Nocera, 2014). When creating our base map, we used two methods to generate two distinct portions of the base map (Fig. S1). The first portion was restricted to the geographic range of our samples (Fig. 1; ∼92% of range) and was determined using the first principle component value determined solely from the fur concentrations of the 14 elements from our validation data set based on 111 samples. The other 15 samples were used for independent verification of the method (below). Wunder & Norris (2019) characterized this type of model as one that directly links tissue chemistry and geography. Fur principle component values were taken from the sPCA analyses, with values interpolated between the locations of our bats using inverse distance weighting (IDW). However, because this approach only uses fur concentrations from samples collected, we are only able project to the limit of our points on a map which does not encompass the entire range of the species (Fig. 1).

To extend our prediction coverage to the second portion (∼8% of range; Fig. 1) of our base map for which no bat samples were available we related soil element concentration to the fur concentrations within the first portion that covered all sampled points. Using those relationships, we transformed soil concentration outside of the sampling extent to hypothetical fur concentrations (Fig. S1). Using the sPCA loadings determined previously, we then transformed the hypothetical fur to PC1 values. Wunder & Norris (2019) describe this approach as an indirectly linked tissue chemistry and geography method. This method is less accurate as more noise is present in the samples due to elemental concentrations being filtered through the food web from soil to fur. As the two maps (directly linked and indirectly linked tissue and geography) do not overlap at any point, they were then merged to create one continuous joint projection across the US and Canada. There was a slight misalignment between the two portions which was left in place to preserve the more accurate direct fur-to-fur comparison present in the base maps generated. Additionally, we tested all metrics for assignment (below) using both the merged maps and the direct geographic relationship (portion 1) separately to ensure adding the second portion did not impact our results.

When generating our base maps we used a leave-three-out (n-3) approach for cross validation. For each process in creating the base map, n-3 samples were used, leaving 3 samples to be assigned back to the base map and testing the performance of our model. This was repeated 100 times, allowing us to assess how our analysis will translate to an independent dataset. When assigning individual bats back to their origin, we used methods originally used in assigning migratory species using hydrogen isotopes values present in bird feathers and bat fur (Fraser et al., 2012; Tonra, Both & Marra, 2014; Fraser, Brooks & Longstaffe, 2017). We did this using the normal probability density function as given in the following equation (Tonra, Both & Marra, 2014): ${\displaystyle f\left(y^{\ast }|{\mu }_i,\sigma \right)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{\frac{-{\left(y^{\ast }-{\mu }_i\right)}^2}{2{\sigma }^2}}}$

where f(y*—μ_i, σ) is the likelihood that an individual with a given principle component fur value (y*) originated from cell i, μ_i is the predicted sPCA value of cell i, and σ² is the variance (σ equaling the standard deviation) within a sampling area which was set at the standard deviation of residuals from the regression between predicted and actual values (σ² = 9.696). The result of this equation is the probability, given a measured fur PC value, came from the cell in question, ranging from 0 to 1. We then iteratively analyzed every cell (0.16^o × 0.16^o; ∼315 km²) of each base map for all individuals and determined the probability a given individual originated from all cells present on the base map.

To evaluate the quality of our model we estimated the probability of assignment for the location where the sample had been collected and considered a prediction ≥50% as correct. We chose 50% as our cut off because we prioritized accuracy (i.e., a correct prediction) above precision (i.e., reducing area of prediction). However, we also evaluated cut-offs of 66% and 75% accuracy as have been used previously (Tonra, Both & Marra, 2014; Seifert et al., 2018). Finally, due to the volant nature of these species, if the prediction value at the sampling location was below 50%, we checked if a point within 10 km (potential nightly travel distance from roost; Hayes & Wiles, 2013) was above our threshold. If a point within 10km was above our 50% threshold, we also considered this a correct assignment to account for the possibility of nightly bat movement around the capture location.

Independent fur sample validation

For further validation we used a subset of random samples from the original 126 (n = 15) that had not been used in base map generation or the leave-three-out analysis. Our logic was that if samples that are independent from the rest of the analyses can also be accurately sourced to their inferred location of origin, it supports the claim that this method can be applied to all L. borealis and can be used in future analyses. Each sample’s fur trace element concentrations were transformed given the sPCA component values previously determined and assigned to a geographic area following the same protocol above.

Trace element concentrations

Concentrations of trace elements in the fur of eastern red bats are similar to values previously reported in fur of other bats (Table 1; for detailed summary see Table S2; Hickey et al., 2001; Zocche et al., 2010; Flache et al., 2015a; Flache et al., 2015b); in addition, levels of variation (reported as SD of mean concentrations) were also similar. Overall, there is a high level of variation in trace element concentrations present among individual samples with all elements having a standard deviation greater than 50% of the mean (Table 1). This suggests that values for individuals collected from different locations differ greatly, which increases the likelihood of correctly assigning individual bats to source areas.

Linking fur and soil element values

All individual elements showed only weak correlations between element values in fur and in the soil with all R² values < 0.16 and only three showing a statistically significant relationship (P < 0.05) (Table 1). This means that individual elements are unlikely to provide sufficient resolution for sourcing individual bats.

Next, we explored the use of composite measure that used multivariate approaches to combine data from different elements in fur. We used a spatial principle component analysis (sPCA) approach to reduce the number of variables from 14 elements to a set of spatial weighted principle component values. We first used Kaiser criterion to evaluate our principle components (PC) from our sPCA from ‘gwpca’ in R (Gollini et al., 2013), with the recommendation to drop PCs with an eigen value below 1.0. By only retaining PCs with eigenvalues above 1.0, we retain factors that extract at least as much information as an equivalent single variable (Kaiser, 1960). In addition, we examined the percent of variation explained by each principle component value. PC 1 explained 50% of the variation present in our data with PC 2 only adding 12% variation explained (Table 2). PCs 3 and larger also added very little to the cumulative variation explained (Table 2). As a consequence, we chose to only retain PC 1 for further analyses. Examination of the first PC shows Al, Cu, Ba, Mg, and Fe had the largest absolute weights (PCA loadings multiplied by average fur element values) on PC 1 relative to their means with values greater than or equal to 1.

The Monte-Carlo permutation test of the sPCA score was not significant (p = 0.07), likely due to similar elemental concentrations in areas that are geographically separated from each other. This indicates that global spatial structures exhibited by elements are not likely due to random variation. To determine if the PC1 values offer a better way to assign individual bats than single elements we compared predicted to actual PC values using a linear model. The PC1-based model had an R² = 0.5 (p < 0.001) which is much higher than the values for individual elements. Consequently, we used PC1 values as our measurement of element variation in subsequent analyses.

Geographic assignment using trace elements

To assess our ability to conduct the geographic assignment of individual bats, we explored three distinct but related factors: accuracy, precision, and geographic patterns. Accuracy is defined by the ability to correctly assign an individual to their origin. Precision is assessed by the amount of area determined as being probable as an area of origin (less area predicted equals higher precision). Finally, geographic patterns are patterns that follow an ordered arrangement through space (e.g., a north/south gradient).

To explore these aspects of our ability to do assignments for the independent dataset, after generating our leave-three-out assignment, we combined the averaged basemaps into one single basemap surface (Fig. 2). Extending our base maps by relating soil to hypothetical fur to principle component values showed a correlation of 0.24 between direct and indirect base map generation. Further, when creating our base maps we compared the merged raster and direct geographic assignment; both showed similar results (accuracy remained consistent, precision ±1%) and so we used the merged basemap for assignment.

Accuracy

We measured accuracy as the number of individuals with a probability over 0.5 or being within 10 km of an area above 0.5. Other threshold values (0.66 and 0.75) that were tested resulted in poor assignment accuracy (56% and 49% correct) and so were not considered further. The probability for the training dataset to be accurate at a 50% threshold was 81%, and was determined by the number of individuals above 0.5 using a leave-three-out analysis (Table 3). Additionally, 12 of our 15 validation independent samples were either above 50% or within 10 km of a cell above our threshold, yielding an accuracy of 80% (Table 3). Because of bats being able to fly ∼10 km from their roost in a given night to feed, we considered these correct assignments. Therefore, using a threshold of 50%, we correctly predict the sampling location ∼80% of the time. Examples of a correct and incorrect assignment from our cross-validation dataset are shown in Fig. 3 (all assignment maps from cross validation in Fig. S3; same order as Table 3).

Table 3:

Locations and results for 15 independent samples.

Locations and prediction information for 15 test individuals and the training dataset following a leave-three-out approach. We report the PC values for each individual, the prediction probability at the sampling location, and the percent of L. borealis range predicted.

Longitude	Latitude	PC values	Probabilty	Range percent
Training	Dataset	Avg	0.81	0.49
−86.10	37.19	20.43	0.52	0.41
−77.04	38.91	34.10	0.94	0.76
−82.60	38.41	97.99	0.00	0.01b
−96.18	36.10	52.52	1.00	0.21
−84.23	30.91	38.45	0.99	0.72
−76.87	38.78	28.79	0.80	0.71
−82.60	34.92	23.99	0.34a	0.71
−97.04	32.90	15.29	0.74	0.22
−77.14	38.97	26.38	0.51	0.65
−95.23	38.97	42.70	0.95	0.59
−99.11	29.89	50.33	0.95	0.29
−84.44	33.71	92.99	0.00	0.03b
−75.48	41.39	17.40	0.32a	0.44
−96.76	46.88	31.58	1.00	0.75
−90.56	42.59	41.69	0.01	0.62b

DOI: 10.7717/peerj.10082/table-3

Notes:

aDenotes values below 50% threshold, but within 10 km of possible sources (above 50%), and are considered successful.

bDenotes those that were incorrect and therefore not included for deteremining precision.