Limited introgression supports division of giraffe into four species.

All giraffe (Giraffa) were previously assigned to a single species (G. camelopardalis) and nine subspecies. However, multi-locus analyses of all subspecies have shown that there are four genetically distinct clades and suggest four giraffe species. This conclusion might not be fully accepted due to limited data and lack of explicit gene flow analyses. Here, we present an extended study based on 21 independent nuclear loci from 137 individuals. Explicit gene flow analyses identify less than one migrant per generation, including between the closely related northern and reticulated giraffe. Thus, gene flow analyses and population genetics of the extended dataset confirm four genetically distinct giraffe clades and support four independent giraffe species. The new findings support a revision of the IUCN classification of giraffe taxonomy. Three of the four species are threatened with extinction, and mostly occurring in politically unstable regions, and as such, require the highest conservation support possible.

INTRODUCTION

Traditionally, giraffe, with their long necks and striking silhouettes (Figure 1), were classified as a single species (Giraffa camelopardalis) with up to eleven subspecies proposed (Lydekker, 1904). However, until recently, the classification into nine subspecies was generally the most accepted one (Dagg & Foster, 1976). It has been shown that in captivity some giraffe subspecies hybridize (Gray, 1972; Lackey, 2011; Lönnig, 2011), which supported the traditional single species concept for giraffe. However, multi‐locus analyses of wild giraffe nuclear loci identified four monophyletic, distinct and evolutionary old groups of giraffe that should be recognized as four distinct species (Fennessy et al., 2016). This finding conflicts with former classifications and has been questioned based on the limited interpretation of traditional data for example, pelage pattern, number of ossicones and geographic distribution (Bercovitch et al., 2017). The initial findings of four giraffe species (Fennessy et al., 2016) could, however, be criticized because it did not involve explicit gene flow analyses. Gene flow analyses are imperative to understanding speciation from a genetic perspective, especially since reproductive isolation is the keystone of the biological species concept (BSC), which is one of the most widely applied (Coyne & Orr, 2004; Mayr, 1942). The BSC implies that there is no or only very limited gene flow between species, but currently, it lacks a clear definition how to deal with gene flow. It has been proposed that one or a limited number of effective migrants (up to 10) per generation (N e m) avoids genetic differentiation of populations and escapes a substantial loss of genetic diversity for neutral traits (Lacy, 1987; Mills & Allendorf, 1996; Vucetich & Waite, 2000; Wright, 1969). Thus, it may be a conservative estimate that limited gene flow of <1 migrant per generation (N e m < 1) can lead to speciation, despite the occurrence of hybridization. It should be noted, that the number of effective migrants is not the same as the actual number of migrating individuals, but an abstract value correlated to the effective population size N e. The N e is typically much smaller (~10%) than the census population size and is equivalent to the size of an idealized population that loses genetic diversity at the same rate as the real population (Frankham, Ballou, & Briscoe, 2010).

Figure 1

Reticulated giraffe (Giraffa reticulata) in the Samburu National Reserve, Kenya (©GCF)

The BSC definition might need to be amended as some species naturally hybridize in the wild and produce fertile offspring for example, bears (Arnold, 2016; Kelly, Whiteley, & Tallmon, 2010; Kumar et al., 2017) and whales (Bérubé & Aguilar, 1998; Spilliaert et al., 1991), and divergence can occur under genetic exchange (Arnold, 2016). While the distinction of four giraffe species is consistent with population genetic analyses (Fennessy et al., 2016), gene flow among giraffe species has not yet been sufficiently analyzed. Here, we revisit the hypothesis of four giraffe species using population genetic methods that explicitly involve gene flow analyses with an increased dataset of 21 nuclear loci and 137 giraffe individuals from 21 locations across Africa (Figure 2; Supporting information Table S1).

Figure 2

Map of Sub‐Saharan Africa with giraffe (sub)species distributions and sampling locations. Geographic ranges (colored shadings) of giraffe as identified by the Giraffe Conservation Foundation (2017) were plotted on a map of Sub‐Saharan Africa. Numbered circles represent sampling locations (for details see Supporting information Table S1). Species and common names as per Fennessy et al. (2016)

MATERIALS AND METHODS

Sampling and DNA Extraction

Tissue samples from all giraffe species and five subspecies were collected by the Giraffe Conservation Foundation (GCF) and partners using remote biopsy darts with country‐specific research permits between 2009 and 2016 in accordance with ethical guidelines and regulations of the respective governments and institutions. All samples were stored either in RNAlater (Invitrogen) or in >95% ethanol. Additional sequences of giraffe individuals were added to the mtDNA dataset of Fennessy et al. (2016) resulting in a total number of 217 giraffe. The geographical origins and individual IDs are shown in Supporting information Table S1. Sample locations and geographical distributions are shown in Figure 2. Additional southern giraffe individuals were included only if they were from a hitherto unrepresented region. DNA was extracted using either a Macherey‐Nagel NucleoSpin Tissue Kit or a standard phenol‐chloroform extraction method. All experimental protocols are in compliance with the guidelines for the best ethical and experimental practices of the Senckenberg Society, as well as with national guidelines of the respective countries.

Amplification and sequencing

We PCR amplified and sequenced the seven intron sequences previously published (Fennessy et al., 2016) for 32 new individuals and developed 14 additional intron sequences as described (Fennessy et al., 2016). The 14 new intron sequences were amplified and sequenced for a total number of 137 individual giraffe and the okapi (Okapia johnstoni). The putatively independent and neutral 21 nuclear gene loci are on different chromosomes or are widely separated from each other in the bovine genome, a close relative with available chromosome level genome data (Supporting information Table S2). PCRs were performed with 10 ng genomic DNA and giraffe and okapi specific primers. Supporting information Table S2 details primer sequences and PCR conditions. We also amplified and sequenced the mitochondrial cytochrome b and control region for all new individuals as described previously (Bock et al., 2014). Each PCR was examined using agarose gel electrophoresis on a 1% agarose gel with ethidium bromide.

Sanger sequencing was performed for the forward and reverse strand using the BigDye terminator sequencing kit 3.1 (Applied Biosystems) with 5 ng of PCR product for each reaction and analyzed on an ABI 3730 DNA Analyzer.

The sequences were manually edited and aligned in Geneious v.6.1.8 (Kearse et al., 2012). Heterozygous insertions/deletions of nuclear sequences were resolved manually or using Indelligent v.1.2 (Dmitriev & Rakitov, 2008) and verified by allele‐specific primers, as necessary. PHASE implemented in DnaSP v.5.10.01 (Librado & Rozas, 2009) was used to derive the allele haplotypes of the nuclear sequences using a threshold of 0.6 and allowing for recombination. All analyses, except of a mtDNA tree analysis, were performed using phased nuclear allele sequences.

Tree analyses

The mitochondrial cytochrome b and control region sequences of 217 giraffe, including new and published sequences (Bock et al., 2014; Brown et al., 2007; Fennessy et al., 2016; 2013; Hassanin et al., 2012; 2007; Winter, Fennessy, Fennessy, & Janke, 2018; Supporting information Table S1) were aligned and concatenated, and a Bayesian analysis was conducted in BEAST v.2.4.5. (Bouckaert et al., 2014). We used the HKY model of sequence evolution (Hasegawa, Kishino, & Yano, 1985), as suggested by jModelTest v.2.1.1 (Darriba, Taboada, Doallo, & Posada, 2012), a log‐normal relaxed clock with 109 generations and sampled every 20,000th iteration. Sequences of two okapis were used as an outgroup.

A multi‐locus Bayesian phylogenetic tree of the 21 intron sequences for 137 individuals and the okapi as outgroup was generated with the StarBEAST2 (Ogilvie, Bouckaert, & Drummond, 2017) package in BEAST v.2.4.5. (Bouckaert et al., 2014, p. 2) under the JC model of nucleotide evolution suggested as best fitting model by jModelTest v.2.1.1 (Darriba et al., 2012). A log‐normal relaxed clock was used with 109 generations and sampling every 20,000th iteration. Convergence of the MCMC runs was analyzed with Tracer v.1.6.0 (Rambaut, Suchard, Xie, & Drummond, 2014), and TreeAnnotator v.2.4.5 (Rambaut & Drummond, 2016) was used to construct a maximum clade credibility tree with 30% burn‐in for both the mtDNA and the nuclear sequences.

Population genetic analyses

Haplotype information for each locus deduced by DnaSP (Librado & Rozas, 2009) was used to code each individual. The haplotype matrix was then used to infer admixture using the Bayesian clustering algorithm implemented in STRUCTURE v.2.3.4. (Pritchard, Stephens, & Donnelly, 2000). For the maximum number of populations (K) between 1 and 10, we sampled 250,000 steps following a 100,000‐step burn‐in, with 40 replicates each. The CLUMPAK webserver (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015) was used to average the results and to infer the most likely K based on the posterior probability of K (Pritchard et al., 2000) and ΔK (Evanno, Regnaut, & Goudet, 2005). Additionally, the most likely K was deduced by eye based on the plot of estimated Ln probability of data (Ln Pr(X|K)) for K between 1 and 10 as described in Pritchard, Wen, and Falush (2010), generated using Structure Harvester (Earl & vonHoldt, 2012). Principal Component Analyses (PCAs) were performed with the R package adegenet (Jombart, 2008) in R v.3.2.3 (R Core Team, 2015) to assess the degree of similarity between defined population scenarios. Pairwise fixation index (F ST) values were calculated in Arlequin v.3.5.2.1 (Excoffier & Lischer, 2010) based on the nuclear haplotypes.

Gene flow analyses

Long‐term average gene flow among and within the giraffe species was calculated in the coalescent genealogy sampler MIGRATE‐N v.3.6.11 (Beerli, 2006; Beerli & Felsenstein, 2001) by estimating the mutation‐scaled population sizes (Θ) for each population and migration rates (M) for each direction between a pair of populations. We used the Brownian motion mutation model and the Bayesian inference analysis strategy, as some parameter combinations are better estimated using the Bayesian approach compared to the Maximum‐likelihood approach (Beerli, 2012). The transition/transversion ratio was set to 2.31 as estimated in MEGA v.7.0.16 (Kumar, Stecher, & Tamura, 2016) based on a concatenated alignment of all 21 loci. Variable mutation rates were considered among loci. We used the default settings for the Θ uniform priors and adjusted the M uniform priors (0; 5,000; 10,000; 1,000) because the upper prior boundary appeared to be too small in initial analyses. Several short runs were performed to check for convergence of the runs. Three long‐chain runs were performed for six million Markov Chain Monte Carlo (MCMC) iterations (60,000 recorded steps) and a burn‐in of 600,000 iterations. An adaptive heating scheme was used with four chains and temperatures set by default with a swapping interval of one. Convergence of the runs was evaluated by the posterior distributions, Effective Sample Size (ESS) with a threshold of ESS ≥ 5,000, and consistency of results between runs. In addition, we estimated short‐term gene flow, as well as the probability of recent hybridization for each individual in BayesAss v.3.0.4 (Wilson & Rannala, 2003) using 100 million MCMC iterations, a burn‐in of 10 million and a sampling interval of 1,000 iterations. Mixing of the chain was improved by adjusting the acceptance rates for proposed changes to the parameters (allele frequencies and inbreeding coefficient) by adapting the mixing parameters for allele frequencies (ΔA) and inbreeding coefficients (ΔF) to 0.30. Convergence was checked in Tracer v.1.6.0 (Rambaut et al., 2014) and by consistency of results of several runs with different initial seeds. Results for short‐term gene flow were visualized in circos plots using the Circos Table Viewer v.0.63‐9 (Krzywinski et al., 2009).

Calculation of gene flow rate

We calculated the N e m using the coalescent‐based estimates for the mutation‐scaled population size Θ and the mutation‐scaled immigration rate M derived from MIGRATE‐N. For autosomal markers, Equation (1) expresses the relationship between Θ (population size of the population receiving migrants) and M ij (corresponding migration rate; Marko & Hart, 2011):

RESULTS

A Bayesian multi‐locus tree analysis of the 21 nuclear loci (total of 16,969 nucleotides) for 137 giraffe, including all traditionally recognized giraffe subspecies (Figure 3a), implies a clear separation into four giraffe clades: (a) a northern giraffe cluster including West African (G. c. peralta), Kordofan (G. c. antiquorum), and Nubian giraffe (G. c. camelopardalis) which includes the former Rothschild′s giraffe (G. c. rothschildi), (b) the reticulated giraffe (G. reticulata), (c) the Masai giraffe (G. tippelskirchi) including former Thornicroft′s giraffe (G. c. thornicrofti), and (d) a southern giraffe cluster (G. giraffa) including Angolan giraffe (G. g. angolensis) and South African (G. g. giraffa). The monophyly of each of these four clades is supported by a posterior probability of p ≥ 0.95. However, in the analyses, the exact relationships of southern and Masai giraffe relative to northern and reticulated giraffe could not be determined with significant probability (p ≈ 0.81) (not shown).

Figure 3

Nuclear phylogeny and population structuring of giraffe. (a) Bayesian multi‐locus tree from 21 nuclear loci and 137 giraffe individuals reconstruct four significant supported (p ≥ 0.95) giraffe clades, corresponding to the four giraffe species (Fennessy et al., 2016). The okapi is used as the outgroup. The asterisks indicate branches with statistical significant support (p ≥ 0.95). The red frame indicates the potential hybrids. (b) STRUCTURE analysis of the dataset, excluding the okapi. The colors indicate the membership in a cluster for each sampling location and individual. K = 4 shows four well‐resolved groups and is supported as best fitting number of clusters by several statistical methods (see Supporting information Figure S2). The grouping into four clusters is consistent with the Bayesian multi‐locus analysis: yellow: northern giraffe, orange: reticulated giraffe, green: Masai giraffe, and blue: southern giraffe. Three individuals within the reticulated giraffe cluster (red arrowheads) indicate potential hybridization with admixture from the northern giraffe. K = 3 merges northern and reticulated giraffe, and at K ≥ 5 no further clustering is evident. (c) PCA axes 1‐2 and axes 1‐3 for four distinct giraffe clusters (1: northern; 2: reticulated; 3: Masai; 4: southern). Colors as in Figure 3b. The 95% confidence intervals are shown as oval outlines. Note that the nonoverlapping confidential intervals in the PCA axes 1‐2, as well as, axes 1‐3 indicate significantly different clusters. Potential hybrids are indicated by black circles. Note. The drawing by Jon B. Hlidberg shows a Nubian giraffe

A mtDNA Bayesian tree (Supporting information Figure S1) confirms the reciprocal monophyly of six distinct subspecies clusters with posterior probability ≥0.95 (Bock et al., 2014; Brown et al., 2007; Fennessy et al., 2016). MtDNA does not support two subspecies of Masai giraffe, even though they appear paraphyletic in the tree, because individuals, which are designated Masai giraffe individuals disrupt a possible reciprocal monophyly and the branch dividing the Masai giraffe into two clusters is not significantly supported. For reticulated giraffe, three individuals do not group as expected but rather fall within the northern giraffe, indicating possible hybridization. However, two of these individuals are from zoos, where hybridization can occur, and have an unknown breeding history. The third individual (LWC01) is a wild giraffe from a geographic range adjacent to the northern giraffe and is a possible natural hybrid. During tissue sampling, this individual was identified phenotypically as reticulated giraffe (J. Fennessy pers. obs.).

Multi‐locus population STRUCTURE analyses (Pritchard et al., 2000) of 21 nuclear loci (Figure 3b) proposes the best clustering into four distinct populations (optimal K = 4) based on the graphical display. At K = 3 the analyses merge the reticulated and the northern giraffe and at K ≥ 5 the analyses do not produce further clustering. Three different statistical methods to interpret the STRUCTURE results (Evanno et al., 2005; Pritchard et al., 2000, 2010) confirm K = 4 being the best fitting number of populations (Supporting information Figure S2). These four clusters conform to the four giraffe clades identified by tree analyses. Intriguingly, STRUCTURE also identifies three potential hybrids between the northern and reticulated giraffe within the reticulated giraffe clade (Figure 2a,b). The distinctness of four unique giraffe clades is in addition supported by Principal Component Analyses (PCAs; Figure 3c) with significant nonoverlapping 95% confidence intervals. PCAs using groups of the seven mtDNA clades do not find more than four distinct clusters (Supporting information Figure S3). Finally, pairwise fixation indices (F ST) of ≥0.237 (statistically significant at p < 0.001) are consistent with the four distinct clusters of giraffe in the tree analyses (Supporting information Table S3).

Separate PCAs and STRUCTURE analyses for each species (Supporting information Figures S4–S7) indicate population substructure within northern giraffe, and potentially in the Masai giraffe, but no further population substructure in southern and reticulated giraffe. Within northern giraffe, STRUCTURE analyses and PCAs up to four clusters can be identified. However, the sample sizes of some populations (three for Ethiopia) are arguably insufficient to draw definitive conclusions if there are further clusters. Within the Masai giraffe, STRUCTURE and PCAs identify potentially two separate clusters, indicating a possible separation of the two geographically most distant populations that have been analyzed for nuclear intron sequences to date. Consistent with the STRUCTURE and PCA analyses, pairwise F ST analyses within each giraffe species find a high level of population differentiation within northern and possibly Masai giraffe, and little differentiation within southern and reticulated giraffe (Supporting information Table S4).

We estimated long‐term gene flow between all four giraffe clades, as well as among subspecies within each of the giraffe species, which show population substructure in STRUCTURE and PCAs using MIGRATE‐N (Beerli, 2006; Beerli & Felsenstein, 2001). All parameters had an ESS > 5,000. Assuming similar mutation rates among all giraffe species, the mutation‐scaled population size theta (Θ) estimates for the four species suggest that the effective population size (N e) is smaller in southern giraffe and Masai giraffe than in northern and reticulated giraffe (Supporting information Table S5a). Thus, the population size in the northern and reticulated giraffe had been larger in the past considering their current numbers (Giraffe Conservation Foundation, 2017). The calculated effective numbers of migrants per generation or gene flow rate (N e m) based on Θ and the mutation‐scaled migration rate (M) (Supporting information Table S5b) indicate generally very low level of gene flow among most of the four giraffe clades with a maximum of one migrant per five generations (N e m ≤ 0.179), with one exception. A higher N e m occurs between the northern and reticulated giraffe with nearly one migrant per generation in the direction of the reticulated giraffe (N e m = 0.945), but much less migration is observed in the opposite direction to northern giraffe (N e m = 0.179). There is little (ca. one in ten) directional gene flow from Masai to reticulated giraffe (N e m = 0.107) and from southern to reticulated giraffe (N e m = 0.104) with nearly zero gene flow in the opposite direction. The gene flow rates for all other species pairs are very low (N e m < 0.065). Within species long‐term gene flow rates are on average higher (N e m > 1) (Supporting information Table S6b). However, between some subspecies, gene flow is also limited, in particular, the geographically isolated West African giraffe (WA).

Finally, short‐term migration rates (m) estimated with BayesAss (Wilson & Rannala, 2003; Figure 4; Supporting information Table S5b) confirm low levels of gene flow among the four giraffe species for the past three generations. The highest migration rates occur from northern, Masai and southern giraffe in the direction of reticulated giraffe, which is expected due to adjacent ranges. The data suggest that approximately 2% (m = 0.021) of the reticulated giraffe population are derived from each of these neighboring species. In comparison with the other gene flow analyses, BayesAss identifies somewhat higher recent migration rates (m) among subspecies within species (Supporting information Table S6b). This is consistent with the lack of genetic differentiation identified by PCA and F ST analyses. The recent migration rates estimated by BayesAss analyses suggest directional gene flow between West African and Kordofan giraffe (m = 0.064), and find gene flow between South African and Angolan giraffe (m = 0.052). Most importantly, however, is that BayesAss does not find any first or second generation hybrids.

Figure 4

Circular migration plot of recent migration rates among four giraffe clades. Recent directional migration rates (m) as estimated by BayesAss and indicated by ribbons connecting one species to another. The color coding of the four species is according to the STRUCTURE clusters (Figure 3b). Peripheral concentric stack bars show relative migration rates in percent. Whereas the inner stack bar shows the outgoing ribbon sizes, the middle stack bar the incoming ribbon sizes and the outer stack bar the combination of both

DISCUSSION

Morphology, ecology, and genetic analyses suggest that there are more than one giraffe species (Brown et al., 2007; Fennessy et al., 2016; Groves & Grubb, 2011; Thomassen, Freedman, Brown, Buermann, & Jacobs, 2013). Here we expand our previous dataset three‐fold and improve the sampling of northern and reticulated giraffe to further study if there is indeed more than one species. The new data allow for the first‐time detailed gene flow and migration analyses with an appropriate amount of data. Among the four giraffe species, gene flow and migration are very limited. As such, the new analyses of the extended nuclear data corroborate the identification of four genetically distinct giraffe species (Fennessy et al., 2016).

Several attempts have been made to define a species, but an unequivocal consensus has not yet been reached (Coyne & Orr, 2004; De Queiroz, 2007). The most commonly applied model is the BSC, which suggests that reproductive isolation is essential to delineate species (Dobzhansky, 1970; Mayr, 1942). By contrast, subspecies or evolutionary significant units are often debated distinctions within a species. Reproductive isolation is also a cornerstone of other species concepts that define species as distinct evolutionary units with limited gene flow to other such units (Avise & Ball, 1990). Therefore, analyzing gene flow among species is a central analysis to delineate species, especially if, like in giraffe, they possibly hybridize in nature. It has been suggested that gene flow among species must be limited to allow genetic differentiation. A value below one migrant per generation (<1 N e m) may be a conservative estimate (Wright, 1969), even if other studies are more liberal and suggest that gene flow rates of <5 N e m (Lacy, 1987) or even <10 N e m (Mills & Allendorf, 1996; Vucetich & Waite, 2000) can allow genetic differentiation and consequently speciation.

The initial finding of four giraffe species was unexpected (Fennessy et al., 2016), as some giraffe interbreed in captivity (Gray, 1972; Lackey, 2011; Lönnig, 2011), and they are highly mobile in the wild (Flanagan, Brown, Fennessy, & Bolger, 2016), both processes which would facilitate admixture. However, we do not observe this in the four giraffe species. If giraffe was in fact one species, mathematical models suggest that long‐term effective gene flow rates in excess of 1–10 migrants per generation would be required to avoid differentiation between populations (Lacy, 1987; Mills & Allendorf, 1996; Vucetich & Waite, 2000; Wright, 1969).

Yet, even among the closely related and neighboring northern and reticulated giraffe lower long‐term effective gene flow rates, <1 N e m, and substantial genetic differentiation are observed, which is consistent with being genetically differentiated species. In addition, among all 137 individuals from a wide geographic distribution, only one natural hybrid was genetically identified. The rare occurrence of hybrids in the wild hints to prezygotic reproduction barriers, because successful hybridization in captivity excludes postzygotic barriers. This further supports the existence of four giraffe species. So far only one possible mechanism for reproductive isolation has been published. Reproduction seems to be synchronized to geographically distinct seasonal rainfall cycles and may contribute to reproductive isolation between the three giraffe species in East Africa (Thomassen et al., 2013).

Population genetic analyses, such as STRUCTURE and PCA of the data set, support the results from the gene flow analyses. The new results are inconsistent with past suggestions of one or possibly six or seven distinct giraffe species (Brown et al., 2007). These previous results were based on nonstringent conclusions from STRUCTURE analyses with 11 separate genetic clusters at K = 13 based on 14 microsatellites, and results of six to seven giraffe clusters based on mtDNA phylogeny (Brown et al., 2007): “11 of the 18 sampling localities resolved as distinct genetic clusters at K = 13,” however, the authors concluded that only “the seven lineages that are reciprocally monophyletic in the mtDNA tree need to be considered evolutionary significant units if not species.” Other findings of up to eight giraffe species were proposed based on a combination of limited genetic analyses (Brown et al., 2007; Hassanin, Ropiquet, Gourmand, Chardonnet, & Rigoulet, 2007) and morphological characteristics (Groves & Grubb, 2011); however, the sample locations of some samples were inaccurate.

Both, Fennessy et al. (2016) and Bock et al. (2014) suggested to subsume Rothschild's giraffe (MF) into the Nubian giraffe, as well as Thornicroft's giraffe (LVNP) into the Masai giraffe because they lack differentiation at mtDNA sequences. Evolutionary differentiation of populations is often first evident in mtDNA because theory suggests that this locus, due to its maternal inheritance and nonrecombining nature, reaches fixation 4‐times more rapidly than nuclear loci (Zink & Barrowclough, 2008). Such population differentiation processes have been reported in natural population of bears (Hailer et al., 2012), humpback whales (Palumbi & Baker, 1994) and macaques (Melnick & Hoelzer, 1992).

While the current mtDNA analyses support previous findings (Fennessy et al., 2016) of Thornicroft's giraffe being subsumed into the Masai giraffe, new and extended nuclear gene datasets identify subtle substructure among them. We emphasize, however, that the nuclear loci have only been sampled from across a limited distribution of the Masai giraffe (Fennessy et al., 2016). Additional sampling of intermediate Masai giraffe populations and additional nuclear gene loci will be necessary to yield more definite results. The first detailed mtDNA analyses on Thornicroft's giraffe (Fennessy, Bock, Tutchings, Brenneman, & Janke, 2013) proposed that while they are not reciprocal monophyletic, the geographic location in Zambia's Luangwa Valley is unique and should, for conservation efforts, tentatively maintain its subspecies status as Thornicroft's giraffe within Masai giraffe (Giraffa tippelskirchi thornicrofti).

Within the northern giraffe, some substructure is evident in PCAs and STRUCTURE analyses for nuclear sequences (Supporting information Figure S6). However, the West African giraffe is a geographically very isolated and small population of ~600 individuals. The geographic distinction between the former Nubian and Kordofan giraffe is unclear and current data suggest that they are not genetically isolated (Fennessy et al., 2016).

With little more than 5,000 northern giraffe, < 15,000 reticulated giraffe and ~34,000 Masai giraffe remaining in the wild (Giraffe Conservation Foundation, 2017), recognizing these—and the southern giraffe—as separate species has an impact on giraffe conservation. Their decline in numbers over the last thirty years (three generations)—northern giraffe (~95%), reticulated giraffe (~60%) and Masai giraffe (~52%)—highlight that these species are threatened with extinction (IUCN, 2017). Giraffe, as a single species, and not four, were recently listed as “Vulnerable” on the IUCN Red List (Muller et al., 2016). The mounting evidence of four giraffe species that now also includes gene flow analyses requires a re‐evaluation of the current IUCN giraffe taxonomy. A higher threat category may be granted to increase conservation management actions.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

AJ, JF, and SW designed and conceived the study. AJ and JF funded the project, and JF collected the samples and provided biological data. SW developed markers, generated and analyzed the data. SW and AJ wrote the manuscript with input from JF.

DATA ACCESSIBILITY STATEMENT

The authors declare that all the data supporting the findings of this study are available within the article and its Supporting Information files. Sequences generated during the study are available at GenBank (https://www.ncbi.nlm.nih.gov/nucleotide/) under the Accession Numbers MG257948–MG262301.

Click here for additional data file.