Vocalizations often vary in structure within a species, from the individual to the population level. Vocal differences among social groups and populations can provide insight into biological processes such as vocal learning and evolutionary divergence, with important conservation implications. As vocal learners of conservation concern, intraspecific vocal variation is of particular interest in elephants. We recorded calls from individuals in multiple wild elephant social groups in two distinct Kenyan populations. We used machine learning to investigate vocal differentiation among individual callers, core groups, bond groups (collections of core groups), and populations. We found clear evidence for vocal distinctiveness at the individual and population level, and evidence for much subtler vocal differences among social groups. Social group membership was a better predictor of call similarity than genetic relatedness, suggesting that subtle vocal differences among social groups may be learned. Vocal divergence among populations and social groups has conservation implications for the effects of social disruption and translocation of elephants.

Data collection

We recorded rumbles from wild adult female elephants in Amboseli National Park, Kenya (“Amboseli”) from 1986–1990 and 1997–2006 and in Samburu and Buffalo Springs National Reserves, Kenya (“Samburu”) from in Nov 2019–Mar 2020 and Jun 2021–Apr 2022. These two populations are 390 km apart with no current gene flow between them due to intervening urban development. Both populations have been continuously monitored for decades, and all individuals can be individually identified by external ear morphology. We focused on adult females (10+ years of age) to ensure that any vocal differences between populations or social groups were not an artifact of age or sex. Our final dataset included calls from 21 adult females in Amboseli (mean ± SD age = 26.2 ± 12.9 years) and 81 adult females in Samburu (mean ± SD age = 25.6 ± 11.2 years). The field recording methods for this dataset have been previously published.

We recorded the identity of the caller and the behavioral context of each call. The caller was identified using behavioral and contextual cues, such as an open mouth, flapping ears, or being the only individual who was not a young calf in the immediate vicinity (calf calls are easily distinguished from adult/subadult calls due to their higher frequency and shorter duration). We only included in the analysis calls for which we were able to identify the caller with certainty. Behavioral context was originally scored using slightly different ethograms in Amboseli and Samburu. To facilitate comparison between these two datasets, we concatenated behavioral context into 9 categories shared across both populations (Supplemental Table S1).

Definition of social groups

We determined the group membership of the elephants in Samburu following a previously published protocol. In brief, we calculated simple ratio association indices between all adult females in the population using observational data collected between January 2019–April 2022, excluding individuals that were seen less than 20 times during this period. We performed a Ward’s hierarchical cluster analysis on association indices, plotted the cumulative number of bifurcations as a function of bifurcation distance, and identified the most significant knot by visual inspection of the plot. We cut the dendrogram at the height corresponding to this knot and designated all clusters below this height as separate core groups. To identify bond groups, we repeated the procedure using only the matriarch (oldest female) of each core group. We did not include Amboseli data in the analysis of core groups and bond groups because we did not have comparable association data for Amboseli, and because 92% of our Amboseli recordings came from individuals belonging to the same (subjectively defined) core group. Data analysis was conducted using R v. 4.2.2.

Call measurement

We measured 94 acoustic features on each call describing the distribution of energy across time and frequency in the mel spectrogram (Supplemental Methods, Supplemental Table S2). A mel spectrogram is similar to a traditional spectrogram (raster plot with time on the x-axis, frequency on the y-axis, and amplitude indicated by pixel darkness) but with frequency transformed to the logarithmic mel scale. While the mel scale was designed to approximate human hearing sensitivity, most other mammals, including elephants, perceive frequency on a similar logarithmic scale.

While the measurements we took from the mel spectrogram capture more of the variation in the calls than traditional measurements such as fundamental frequency and formants, they also capture background noise and thus may be more susceptible to influence from the recording equipment. To ensure that differences between Amboseli and Samburu could not be attributed to the different recording gear used in each population, we also traced the second harmonic (f₁) contour of each call, which is highly robust to recording equipment differences. We calculated nine summary statistics of the f₁ contour (Supplemental Table S2).

Individual, social unit, and population assignment accuracy

To test the hypothesis that Amboseli and Samburu elephants exhibit population-level acoustic differences, we ran a random forest (500 trees, 6 variables/node, 60% of observations/tree, minimum node size = 1, no maximum tree depth) to predict population as a function of the mel spectrogram acoustic features. A random forest is a machine learning model comprised of many decision trees, each using a randomly selected fraction of the data. This analysis included 938 calls from 81 individuals in Samburu, and 414 calls from 21 individuals in Amboseli. As random forests are biased towards more numerous classes, we balanced the dataset by randomly subsampling the Samburu observations so there were an equal number of observations in each class. To ensure that the model could only predict population using cues that generalized across callers, we randomly selected 20% of the callers with at least five calls each from each population and allocated all calls from these callers to the test set, with the remaining calls allocated to the training set. We calculated the proportion of observations in the test set that were classified correctly (classification accuracy) and ran a one-tailed exact binomial test comparing the classification accuracy to the proportion that would have been classified accurately if the model always guessed the most common group in the training set (“majority” or “zero rate” classifier). The majority classifier is a standard baseline model used in machine learning. When there are the same number of observations in each class, the majority classifier is mathematically equivalent to the weighted expectation (guessing each class with a probability equal to the proportion of the data comprised by that class). However, when the data are imbalanced, the majority classifier always outperforms the weighted expectation. We repeated the above process 10,000 times and calculated the median P-value across all runs. The number of calls allocated to the test set varied across runs because different callers produced different numbers of calls. The mean ± SD proportion of the calls allocated to the test set across 10,000 runs was 0.16 ± 0.05. Note that although the full dataset was balanced by subsampling an equal number of calls from each population, the data were not perfectly balanced within the training and test sets, because the number of calls from each population allocated to the test set depended on the number of calls per caller. Moreover, the majority population in the training set was not always the majority population in the test set. Thus, the majority classifier accuracy could be greater or less than 50%.

To determine if vocal differences between the two populations could be an artifact of the different recording equipment used in each population, we ran the same model using the f₁ contour measurements instead of the mel spectrogram measurements. We also ran a logistic regression model with population as the response variable and the f₁ contour measurements and caller age as the regressors, to determine which acoustic features had a significant relationship to population.

To test the hypotheses that elephants exhibit vocal signatures of group identity at the bond group or core group level, we ran two additional random forest models (same hyperparameters) to predict bond group and core group, respectively, as a function of the mel spectrogram acoustic features, using only data from Samburu. To ensure that the bond group model could only use acoustic features that generalized across the entire bond group, rather than features specific to core groups, to predict bond group, we restricted the dataset to bond groups that contained at least two core groups with at least five calls each in our dataset (6 bond groups, 931 calls). For each iteration of the model, we randomly selected one core group from each bond group and allocated all calls from those core groups to the test set. To ensure that the core group model could only use acoustic features that generalized across the entire core group, rather features specific to individual callers, to predict core group, we restricted the dataset to core groups that contained at least two individuals with at least five calls each in our dataset (7 core groups, 794 calls). For each iteration of the model, we randomly selected 20% of the callers with at least five calls each from each core group and allocated all calls from those individuals to the test set. The data was severely imbalanced across bond groups and core groups, but we were unable to balance it by subsampling the more numerous classes because doing so would have reduced the sample size excessively. We ran 10,000 iterations for each model, calculating the classification accuracy and p-value for each run as before. The mean ± SD proportion of the calls allocated to the test set was 0.32 ± 0.17 for the bond group model and 0.23 ± 0.08 for the core group model.

To test the hypothesis that elephant rumbles are individually specific, we ran a fifth random forest (same hyperparameters) to predict individual caller identity as a function of the mel spectrogram acoustic features. As calls produced by the same caller on the same date might exhibit similar features due to temporary circumstances such as caller’s internal state, behavioral context, and ambient conditions, we randomly selected one date for each caller and held out all calls from these caller-dates as the test set. We used callers from both populations for this analysis, but only included callers that produced at least three calls on at least two different dates each (427 calls from 15 callers in Samburu, 218 calls from 9 callers in Amboseli). We calculated the classification accuracy and p-value for each of 10,000 iterations as before. The mean ± SD proportion of calls allocated to the test set was 0.20 ± 0.02.

Assessment of call similarity

Elephant rumbles vary with the behavioral context and the individual identity and age of the caller. To determine if there were acoustic differences among populations, bond groups, or core groups that could not be explained by behavioral context, caller identity, or caller age, we calculated random forest proximity scores between each possible pair of calls. The random forest proximity score for a given pair of calls was the proportion of trees for which both calls were classified in the same terminal node, adjusted for the size of the node, and represented a metric of call similarity in terms of the acoustic features most relevant to predicting the response variable. We calculated proximity scores from four different random forests (population with mel spectrogram features, population with f₁ contour features, bond group, and core group), using the same hyperparameters and subsets of the data as before except that we increased the number of trees to 8000 and did not balance the data by subsampling or hold out any observations as a test set (population models: all observations, n=1352; bond group model: all bond groups in Samburu containing at least 2 core groups with at least 5 calls each, n=931; core group model: all core groups in Samburu containing at least 2 individuals with at least 5 calls each, n=794).

For each set of proximity scores, we ran a generalized linear mixed model with a gamma error distribution and a log link function. In each model, pairwise call proximity score was the response variable and pair ID (unique identifier for a given pair of callers) was a random effect. Fixed effects were “pair class” (whether the two calls in a given pair came from the same population, bond group, or core group, depending on the model) and the scaled and centered absolute value of the age difference between the two callers. For the population models, the “pair class” variable was binary: a pair of calls could either be from the same population or not. For the bond group and core group models, the “pair class” variable had three possible states: same core group, different core groups within the same bond group, or different bond groups within the Samburu population. To ensure that differences between populations or social groups were not an artifact of individual identity or behavioral context, we only included pairs of calls with the same behavioral context and different callers. We also only included calls for which we were certain of the behavioral context (see Supplemental Methods). This resulted in a sample size of 149,640 call pairs for the two population models, 97,586 call pairs for the bond group model, and 72,437 call pairs for the core group model. As proximity scores could be 0, we added 0.00001 to all proximity scores so all the values would be positive. For the bond group and core group models, we used the R package “emmeans” to examine the pairwise contrasts for each level of the “pair class” variable.

To assess whether vocal similarity between individuals was better explained by social affiliation or genetic relatedness, we ran two additional gamma regressions on call pairs, modeling call proximity score as a function of “binary pair class” (see below), caller age difference, and caller genetic relatedness. One gamma model used proximity scores extracted from the random forest trained to predict core groups, and the other used proximity scores extracted from the random forest trained to predict bond groups. These two sets of proximity scores represented the pairwise similarity between calls in terms of the features most relevant to predicting core group membership or bond group membership, respectively. For the model using proximity scores extracted from the bond group random forest, binary pair class indicated whether the two calls in a pair were from the same bond group. For the model using proximity scores extracted from the core group random forest, binary pair class indicated whether the two calls in a pair were from the same core group. These gamma models could only be run on the subset of callers for which we had genetic relatedness data (13 individuals for the bond group model, 8 for the core group model). Genetic relatedness was calculated in a previously published study using 20 microsatellite loci extracted from tissue or dung samples collected before 2006. Due to social disruption caused by poaching, many core groups and bond groups in the Samburu population now include unrelated individuals, which made it possible to independently assess the effects of social affiliation and relatedness on call similarity.

All statistical analyses were performed in R v. 4.2.2, and 0.05 was used as the significance threshold for all tests.