Genetic structuring among colonies of a pantropical seabird: Implication for subspecies validation and conservation

Abstract Investigations of the genetic structure of populations over the entire range of a species yield valuable information about connectivity among populations. Seabirds are an intriguing taxon in this regard because they move extensively when not breeding, facilitating intermixing of populations, but breed consistently on the same isolated islands, restricting gene flow among populations. The degree of genetic structuring of populations varies extensively among seabird species but they have been understudied in their tropical ranges. Here, we address this across a broad spatial scale by using microsatellite and mitochondrial data to explore the population connectivity of 13 breeding populations representing the six subspecies of the white‐tailed tropicbird (Phaethon lepturus) in the Atlantic, Indian, and Pacific Oceans. Our primary aim was to identify appropriate conservation units for this little known species. Three morphometric characters were also examined in the subspecies. We found a clear pattern of population structuring with four genetic groups. The most ancient and the most isolated group was in the northwestern Atlantic Ocean. The South Atlantic populations and South Mozambique Channel population on Europa were genetically isolated and may have had a common ancestor. Birds from the Indo‐Pacific region showed unclear and weak genetic differentiation. This structuring was most well defined from nuclear and mtDNA markers but was less well resolved by morphological data. The validity of classifying white‐tailed tropicbirds into six distinct subspecies is discussed in light of our new findings. From a conservation standpoint our results highlight that the three most threatened conservation units for this species are the two subspecies of the tropical North and South Atlantic Oceans and that of Europa Island in the Indian Ocean.


| INTRODUC TI ON
Widely distributed species consist of distinct populations that are variously connected to each other though current and historic overlaps in ranges, resulting in different patterns of genetic structuring at global and local scales (Avise & Ball, 1990;Frankham, Ballou, & Briscoe, 2004). In mobile animals, life-history traits such as dispersal and philopatry, as well as geographic isolation, shape population genetic connectivity (Dobzhansky & Dobzhansky, 1970;Greenwood, 1980;Matthiopoulos, Harwood, & Thomas, 2005;Taylor & Friesen, 2012;Wright, 1943). Many species with wide distribution are subdivided into subspecies (Ball & Avise, 1992). The distinction between a subspecies and a species is at time far from clear but a lineage can be considered as a new species when it acquires sufficiently different properties from others. Such differences can include being phenetically distinguishable, diagnosable, reciprocally monophyletic, reproductively incompatible, and ecologically distinct (de Queiroz, 2007).
There is a critical (and ongoing) need for delineation of species and subspecies because it is fundamental in elucidating key processes in ecology (e.g., population connectivity) and conservation biology (e.g., definition of conservation priority units; Dayrat, 2005;Gutiérrez & Helgen, 2013). Traditional morphology-based taxonomy defining "morphotypes" (i.e., taxa described solely from morphological traits) faces challenges when applied to taxa that are either cryptic (i.e., a group of individuals that is morphologically indistinguishable but incapable of interbreeding) or display considerable phenotypic plasticity. Therefore, Zink (2004) and Sackett et al. (2014) argued that only taxa defined by the congruence of multiple morphological or molecular characters should be recognized as distinct subspecies.
Of the three species of tropicbird (Phaethontidae), the whitetailed tropicbird Phaethon lepturus is the most common (Lee & Walsh-McGehee, 1998). This pantropical seabird is widely distributed in the Atlantic, Pacific, and Indian Oceans between 30°N and 30°S (del Hoyo, Elliott, Sargatal, & Christie, 1992). Although its global conservation status is of "Least Concern," the white-tailed tropicbird is decreasing globally predominantly because of predation by invasive species (BirdLife International, 2018). Declining populations are reported in the Atlantic islands where tropicbirds have been depredated by introduced brown rats Rattus norvegicus, black rats R. rattus, feral domestic cats Felis sylvestris catus, dogs Canis lupus familiaris, American Crows Corvus brachyrhynchos (on Bermuda), and tegu 13 breeding populations representing the six subspecies of the white-tailed tropicbird (Phaethon lepturus) in the Atlantic, Indian, and Pacific Oceans. Our primary aim was to identify appropriate conservation units for this little known species. Three morphometric characters were also examined in the subspecies. We found a clear pattern of population structuring with four genetic groups. The most ancient and the most isolated group was in the northwestern Atlantic Ocean. The South Atlantic populations and South Mozambique Channel population on Europa were genetically isolated and may have had a common ancestor. Birds from the Indo-Pacific region showed unclear and weak genetic differentiation. This structuring was most well defined from nuclear and mtDNA markers but was less well resolved by morphological data. The validity of classifying white-tailed tropicbirds into six distinct subspecies is discussed in light of our new findings. From a conservation standpoint our results highlight that the three most threatened conservation units for this species are the two subspecies of the tropical North and South Atlantic Oceans and that of Europa Island in the Indian Ocean.

K E Y W O R D S
conservation status, genetic structure, Phaethon lepturus, subspecies status lizards Salvator merianae (in the Fernando de Noronha archipelago) and where seabirds generally have long endured persecution from humans (Efe, Serafini, & Nunes, 2018;Lee & Walsh-McGehee, 2000;Nunes, Efe, Freitas, & Bugoni, 2017). In Bermuda, one of the main threats was habitat destruction through storm damage; in 2003, a category 3 hurricane damaged 50% of natural nest sites of whitetailed tropicbirds in one of national parks (PT unpublished data). In the Indian Ocean, the breeding tropicbird population of Aride Island, Seychelles, suffered a marked decline of 60% between 1989 and 1998 (Bowler, Betts, Bullock, & Ramos, 2002) which continues to the present day as a 5.4% annual rate of decline (Catry et al., 2009).
Threats to tropicbirds on Aride have not been clearly identified to date but 23.2% of recorded adult mortalities are due to entanglement in Pisonia sticky seeds with birds losing the ability to fly off.
Based on morphological differences, six subspecies of P. lepturus are currently recognized, three occurring in the Indian Ocean, two in the Atlantic Ocean, and one in the Pacific Ocean (del Hoyo et al., 1992;Le Corre & Jouventin, 1999). Three large subspecies breed in the western (P. l. lepturus) and eastern (P. l. fulvus) Indian Ocean and in the northwestern Atlantic Ocean (P. l. catesbyi). The three small subspecies breed in the Pacific (P. l. dorotheae), Indian (South Mozambique Channel; P. l. europae), and south Atlantic Oceans (P. l. ascensionis). The species has golden or apricot plumage morphs in some populations in which birds have white parts of the body plumage that are tinged with these shades. Populations of P. l. fulvus and P. l. europae comprise large fractions of apricot and golden birds, respectively, whereas other subspecies contain mainly white individuals (Le Corre & Jouventin, 1999; see Table 1 and Figure 1).
Such body size and plumage differences among subspecies and locations suggest some level of genetic isolation among populations (Le Corre & Jouventin, 1999). In addition, even neighboring populations could have sharp genetic structuring (Nunes et al., 2017;Wiley et al., 2012). Therefore, the aim of the present study was to explore genetic diversity and population genetic structuring across the entire range of the white-tailed tropicbird and to infer genetic relationships among populations. Both microsatellite and mtDNA markers were used to investigate genetic structure of 13 breeding populations of the Atlantic, Indian, and Pacific Oceans, including birds from all of the six recognized subspecies. We also performed morphometric analyses to examine patterns of phenotypic differentiation among 11 populations of birds for which data were available.
Our findings allow us to discuss subspecies validity and appropriate units for effective conservation strategies.

| Study sites and sampling
Thirteen populations were sampled between 2008 and 2013 in Atlantic, Indian, and Pacific Oceans (see Table 2 and Figure 1 for further details). In each population (except that at Réunion Island), all birds were sampled from the same local area and thus can be considered to belong to the same breeding colony. On Réunion Island, we sampled birds from the local Wildlife Rescue Center (SEOR) that originated from various locations across the island. Details of field researchers and local licenses (when required) are provided in supporting information (Appendix S1, Table S1). Blood was collected from 382 birds by venepuncture of the brachial or metatarsal veins and samples were stored in 70% ethanol or by sterile syringe and needle and stored dry on filter paper.

| Biometrics
Morphometric measurements were obtained from breeding adults in 11 populations (Table 2): (a) wing length (i.e., flattened wing from the carpal joint to the tip of the longest primary), (b) tarsus length (i.e., maximum tarsus length), and (c) bill length (i.e., exposed culmen) following the British Trust for Ornithology's (BTO's) Ringers' Manual (Redfern & Clark, 2001). Morphometric measurements of 616 individuals of white-tailed tropicbirds are available in supporting information (Appendix S1, Table S2). Statistical differences in morphometric measurements of these birds between pairs of focal populations were examined using pairwise Wilcoxon rank sum tests corrected after Benjamini and Hochberg (1995), in R 3.4.2 (R Core Team, 2017). Global patterns in morphometric differences among birds and the populations to which they belong were investigated using a principal components analysis (PCA) performed on a correlation matrix of three morphometric measurements from 421 birds of 10 populations individuals for which we had measurements for all F I G U R E 1 (a) Locations (stars) and associated subspecies of white-tailed tropicbird populations sampled in the Atlantic, Indian and Pacific Oceans (see Tables 1 and 2  three characters (Appendix S1, Table S2) using the ade4 package in R (Dray & Dufour, 2007

| Microsatellite genotyping
Total DNA was extracted from blood or tissue samples using the QIAmp Blood and Tissue kit (Qiagen). Genotyping was conducted for 10 polymorphic microsatellite loci (described in Humeau et al., 2011) on DNA extracts from 382 individuals ( white-tailed tropicbirds are available in the supporting information (Appendix S1, Table S3).

| Mitochondrial DNA sequencing
Mitochondrial DNA variation was analyzed in a subset of samples used for microsatellite genotyping (n = 123, Table 2) and was assayed following the amplification of part of the cytochrome oxidase subunit I (COI) region. Sequences from white-tailed tropicbirds (Kennedy & Spencer, 2004; GenBank accessions AY369055, JN801349) were used to design two primers: CO1Fcons 5′-GACCGAAACCTAAAYACCACA-3′ and CO1Rcons see Appendix S1, Table S4 for details). In addition to the P. lepturus samples, DNA from six red-tailed tropicbirds of P. rubricauda was also included in the phylogenetic analysis as an outgroup (GenBank accession numbers: MT073267 and MT073268; Appendix S1, Table   S4). This species is the phylogenetically closest interspecific relative to P. lepturus (Kennedy & Spencer, 2004).  (Rousset, 2008). The p-values were corrected using the Benjamini and Yekutieli (2001) method for multiple comparisons (Narum, 2006). The mean observed number of alleles per locus (AL) and the number of private alleles per population (AP) were computed using GenAlEx 6.5 (Peakall & Smouse, 2012). Allelic richness (AR-El Mousadik & Petit, 1996), adjusted for discrepancies in sample size by incorporating a rarefaction method as implemented using FSTAT 2.9.3 (Goudet, 2001), was used to make comparisons of the mean number of alleles among populations (except for the populations from Fernando de Noronha and Madagascar that contained only four individuals each). Observed heterozygosity (H O ), unbiased expected heterozygosity estimated according to Nei (1978) (H E ), and Wright's F-statistics F IS according to the method of Weir and Cockerham (1984) were calculated for each population using GENETIX 4.05.2 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 1996-2004. Deviations from Hardy-Weinberg equilibrium (HWE)

| Genetic diversity
were tested for each population using GENETIX 4.05.2.
Heterozygosity excess tests were performed with the program BOTTLENECK 1.2 (Piry, Luikart, & Cornuet, 1999)  Statistical significance of the number of loci with heterozygosity excess as expected in bottlenecked populations (Luikart et al., 1998) was evaluated by a one-tailed Wilcoxon signed-rank test from 10 4 simulation replicates. The L-shaped method illustrates the frequency of rare alleles in the populations; an L-shape graph indicates that the population is in mutation-drift equilibrium (Luikart et al., 1998). The BOTTLENECK software was also used to establish if allele frequency distribution was L-shaped. Finally, test for signatures of a bottleneck by the M-ratio method (Garza & Williamson, 2001) was applied using the StrataG package in R (Archer, Adams, & Schneiders, 2016).
The M-ratio statistic indicates the number of unoccupied potential allelic states and was shown to be small (<0.68 in a dataset with seven loci) when a severe population decline had occurred (Garza & Williamson, 2001).

| Genetic differentiation and structuring
A hierarchical analysis of molecular variance (AMOVA) was performed on microsatellite allele identity and mtDNA to determine how genetic diversity was distributed within and among subspecies and populations. Statistical significance was determined in ARLEQUIN 3.5.1.3 (Excoffier & Lischer, 2010) with 1,000 permutations.
Genetic differentiation among all pairs of populations with sample sizes ≥ 5 was assessed by calculating pairwise F ST and Φ ST values for microsatellite and mtDNA data, respectively. F ST was computed between pairs of populations following Weir and Cockerham (1984) and Whitlock (2011), and statistical significance was tested by 10 4  (Pritchard et al., 2000). For each value (1-13) of the number of independent genetic clusters or K, we ran 10 6 iterations 20 times (after a burn-in of 5 × 10 5 steps). For choosing the optimal number of clusters, three criteria were used; (1) the log likelihood given K (L(K); Pritchard et al., 2000), (2) the second-order rate of change of mean log likelihood (∆K; Evanno, Regnaut, & Goudet, 2005), and (3) the median value of L(K). The first two were calculated using STRUCTURE HARVESTER online Web server (Earl & VonHoldt, 2012). The third was calculated using CLUMPAK software (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015) that was also used to find the optimal individual alignments of replicated cluster analyses and to visualize the results.
Population structure was also explored by performing a  (Jombart et al., 2010). Therefore, the rate of decrease in BIC values was visually examined to identify values of K, after which BIC values decreased only subtly (Jombart et al., 2010). DAPC was applied using the Adegenet package 2.1.1 in R (Jombart, 2008).
A hierarchical analysis of the genetic structuring using Bayesian and DAPC procedures was performed at different spatial scales to detect a possible pattern among the less isolated populations. This was performed: (a) in all populations, (b) after removing the much-differentiated Bermuda population, and (c) after removing the three differentiated genetic groups on Bermuda, Ascension/Fernando de Noronha and Europa.

| Network and phylogenetic analysis
A statistical parsimony haplotype network was obtained from mtDNA sequences (123 individuals sampled in 13 populations) using the pegas package 0.11 in R (Paradis, 2010). The haplotype network was built using an infinite site model (Hamming distance) (Templeton, Crandall, & Sing, 1992) to display the maternal connections among populations and subspecies of P. lepturus.

| Phenotypic variation
Based on three morphometric measurements, birds of the subspecies P. l. catesbyi (i.e., in the Bermuda population) had significantly longer tarsi than birds in the other populations (Table 3). Birds of the subspecies P. l. europae (i.e., in the Europa population) were consistently smaller for all three measured traits compared with the remainder (differences of 25, 6, and 6 mm between Bermuda and Europa populations for wing (all Ps ≤ .006), tarsus (all Ps ≤ 4 × 10 -6 ), and bill (all Ps ≤ 9 × 10 -6 ) lengths, respectively; Table 3). No birds in the other subspecies differed significantly in size in pairwise comparisons (Table 3). The PCA revealed the separation of individuals along the PC1 axis (explaining 59% of the total variance) into three main clusters corresponding to subspecies P. l. catesbyi, P. l. europae, and the others (P. l. ascensionis, P. l. lepturus and P. l. fulvus; Figure 2). The populations were significantly different in morphology when tested by parametric MANOVA (Wilk's lambda = 0.33, df = 9, p < 2.2 × 10 -16 ). Pairwise permutation MANOVAs showed that P. l. catesbyi and P. l. europae were significantly bigger and smaller, respectively, compared to all other populations (all Ps < 0.05).

| Genetic diversity and test of bottleneck
No null alleles, large-allele dropout nor stutter peaks were detected for the 10 microsatellite loci. The percentage of missing data was 0.81%. Linkage disequilibrium among loci was detected for three of the 45 loci pairs (p < .05) but no significant linkage disequilibrium was observed among any of the loci after the Benjamini and TA B L E 3 Mean morphometric measurements (±1 SD) of five subspecies of Phaethon lepturus in 11 populations (see Table 1 for further details)

Subspecies Morphometric measurements
Wing length (mm)

| Genetic differentiation and population structure
The AMOVA using microsatellite allele identity indicated that 65% of microsatellite variation was attributed to differences among the delineated subspecies, 10% among the populations and 25% among individuals (all Ps < 0.001). With mtDNA, 73% of the variance was attributed to differences among subspecies, 13% to among populations and 14% to among individuals (all Ps < 0.001).  Table S6).

| Microsatellite differentiation and clustering
Clustering of microsatellite genotypes using STRUCTURE analysis showed that the best-supported model contained three to eight genetic clusters, depending on the optimal K method used (Figure 4a) and suggested a hierarchical structure. Without F I G U R E 2 Scatterplot from a principal component analysis (PCA) of morphological variation based on three morphometric measurements [wing length (WL), tarsus length (TL) and bill length (BL)] taken from 421 Phaethon lepturus from 10 localities. The total variance explained by PC1 and PC2 was 59% and 25%, respectively. Colors correspond to populations with population codes detailed in Table 2 TL  assignment-success rate) were isolated from the others (Figure 5b).
All the populations of the Indo-Pacific Ocean except Europa were admixed, with a maximum of 71% assignment-success rate between two clusters (Figure 5c).
Phylogenetic analysis confirmed the same relationships as in the haplotype network ( Figure 6). Both phylogenetic analyses

| Correlations among genetic, morphometric, and geographic distances
Genetic distances based on nuclear microsatellites and mtDNA sequences estimated for 10 populations showed significant relation-

| D ISCUSS I ON
The present study has provided data on genetic and morphological differences among putative white-tailed tropicbird subspecies. Patterns of variation in white-tailed tropicbirds revealed significant differences at different geographic scales and sug-

| Morphological variation among subspecies
Morphometric measurements of birds showed clear differentiation among the larger subspecies P. l. catesbyi, the smaller subspecies P.
l. europae and the other subspecies that were intermediate in size (Table 3 and Figure 2). Our results were not entirely consistent with those of Le Corre and Jouventin (1999) who distinguished between large subspecies (i.e., P. l. catesbyi, P. l. lepturus, and P. l. fulvus) and small subspecies (i.e., P. l. ascensionis, P. l. europae, and P. l. dorotheae).
This difference in findings of these two studies may be partly explained by Le Corre and Jouventin (1999) sampling birds of the subspecies ascencionis from the São Tomé population as opposed to the present study that included birds from the Ascension and Fernando de Noronha populations. There are slight morphological differences between birds from the western and eastern South Atlantic, with birds from Ascension/Fernando de Noronha populations slightly bigger than birds from the São Tomé population (see table 2 in Le Corre and Jouventin (1999) and Table 3 in the present study).
As previously reported by Le Corre and Jouventin (1999), body size appeared to decrease from the highest (Bermuda, 32°N) to the lowest (Europa, 22°S) latitudes but variation was only subtle with none of the correlations of any morphological traits with latitude or longitude being statistically significant (data not shown). Body size variation of P. lepturus could be explained, in part, by Bergmann's Rule, applied at an intraspecific level, as suggested by Le Corre and F I G U R E 6 Phylogenetic relationships and estimated divergence times based on COI mtDNA sequences (872 bp) from 123 Phaethon lepturus from 13 populations. Statistical support is given above branches; number shown in order for the respective analysis (Beast posterior probabilities/Maximum Likelihood bootstrap values) only if p > .50 and bootstrap values >50. Numbers on the x-axis refer to millions of years (MYr) since the present; mean estimated times in MYr and the 95% intervals are indicated for each major node corresponding to the four phylogroups. Colors correspond to populations whose codes and associated subspecies names are detailed in Table 2 1.0 0.5  (1999). Phenotypic structuring results from the selection of multiple pressures to which an organism is exposed (Mayr, 1956), such as differences in local oceanographic characteristics around distinct seabird populations (Friesen, 2015). Environmental characteristics (e.g., sea surface temperature, primary productivity) and characteristics relating to foraging behaviour seem to be correlated with the distribution of phenotypes in the marine realm, as has been suggested for seabirds (Weimerskirch, Zimmermann, & Prince, 2001;Yamamoto et al., 2016). Islands are separated by 39° of latitude (or 7,800 km), and so it is likely that latitudinal distance plays a major role in their genetic isolation. There is a pressing need to investigate the genetic structure of birds along the southeast-northwest gradient of the tropical Atlantic (including Great and Lesser Antilles), to identify the spatial limits of northern and southern genetic clusters.

| Genetic structuring among subspecies
The phylogroup of P. l. ascensionis and P. l. europae could have a common ancestor that probably diverged from the ancestor of the Indo-Pacific region about 290 kya ( Figure 6). All nuclear microsatellite analysis showed a strong isolation for each of these two populations. The geographic localization of the ancestor could not be determined and would not have been reliable with only one sequence of <1,000 bps carried a tiny amount of phylogenetic signal. Nevertheless, sporadic gene flow could be possible between Atlantic and Indian Oceans. Some birds of the tropical Atlantic (i.e., from Ascension and Fernando de Noronha, but also from São Tomé and Príncipe in the Gulf of Guinea-not included in our study) may have migrated to the Indian Ocean, probably through the southern tip of Africa. Similar gene flow between Atlantic and Indian Oceans was proposed for petrels (Pterodroma spp.;Booth Jones et al., 2017;Brown et al., 2010) and between Atlantic and Pacific Oceans for magnificent frigatebirds Fregata magnificens (Hailer et al., 2010).  (Bourjea et al., 2007). The red-footed booby Sula sula is a polymorphic pantropical seabird and most populations in the Indian Ocean predominantly consist of the white-tailed white morph except at Europa Island where >95% are of the white-tailed brown morph, suggesting strong isolation (Le Corre & Jouventin, 1999) as confirmed by a recent genetic study (Danckwerts, 2018). We suggest that the strong genetic isolation of P. l. europae on Europa Island may be due to several factors such as natal philopatry and breeding fidelity, breeding phenology, ecological specialization, and sexual selection (Friesen, Burg, & McCoy, 2007;Lombal et al., 2020;Sexton, Hangartner, & Hoffmann, 2014;Uy, Irwin, & Webster, 2018).

Birds of Europa
Surprisingly, the only population of the Pacific Ocean that we included in our study (Tahiti, French Polynesia) was genetically indistinguishable from populations of the Indian Ocean (except birds in the Europa population). The lack of nuclear and mtDNA genetic differentiation between Indo-Pacific subspecies (i.e., P. l. lepturus and P. l. dorotheae, represented by one breeding population in Tahiti) demonstrates that the Indonesian archipelago is not a physical barrier to seabird dispersion, suggesting that gene flow occurs among populations in the two ocean systems. An alternative scenario that would explain this result would be from only recent isolation of a previously panmictic population. An ongoing regional study of migration strategies of 10 seabird species of the western Indian Ocean has shown that some species, such as Lesser frigatebirds (Fregata ariel), sooty terns (Onychoprion fuscatus), and red-tailed tropicbirds (Phaethon rubricauda), occasionally fly over the Timor Sea and the Indonesian straits to reach the Pacific Ocean (Jaeger et al., 2017;Le Corre et al., 2012;Weimerskirch et al., 2017;MLC unpublished data

| Implications for subspecies validity and conservation
Our findings question the validity of classifying white-tailed tropicbirds into six distinct subspecies and help to define appropriate conservation units (Frankham, 2010;Haig et al., 2011;Taylor & Friesen, 2012). The results of AMOVA tests using microsatellite and mtDNA markers confirmed that subspecies designation explained the global genetic variation of the species at large better than did the population designation. Variation among the delineated subspecies was high (i.e., >65% for both markers), but only a low percentage (i.e., 10%-13%) of the variation could be attributed to differences among populations.

| Isolated and threatened subspecies: P. l. catesbyi, P. l. ascensionis, and P. l. europae
Both nuclear and mtDNA data identified three distinct lineages of birds on Bermuda, Ascension/Fernando de Noronha and Europa, corresponding to P. l. catesbyi, P. l. ascensionis, and P. l. europae subspecies, respectively. Each lineage was restricted to its own island with no co-occurrence of lineages on the same island. These three subspecies should be maintained and be defined as a priority for conservation. Low genetic diversity indices (Table 4)  in P. l. europae and P. l. ascensionis at Ascension (Table 4) and at Fernando de Noronha (see Nunes et al., 2017).
The Bermuda genetic group probably includes all populations of the northwest tropical Atlantic Ocean (i.e., P. l. catesbyi). Lee and Walsh-McGehee (2000) reported that the subspecies had declined in regional population size by 50% since 1984, probably due to loss of breeding habitats and predation by invasive mammals. They estimated that the region contained 5,500 pairs, of which 2,500 (or 45%) breed in Bermuda, 1,000 (or 19%) in the Bahamas, 1,000 in Hispaniola and the Greater Antilles and with the remainder (or 17%) breeding in 20 small populations. Because the species is highly philopatric, individuals are unlikely to relocate their breeding sites to avoid further threats (Wingate & Talbot, 2003). Their breeding success (as measured by the percentage of nests regularly visited by adults that fledge a chick) decreased from 66.6% in 1970-1983to 48.8% in 2001(Wingate & Talbot, 2003. The nonsignificant tests of genetic bottleneck could be unreliable due to small sample sizes (Piry et al., 1999). Phaethon lepturus catesbyi should be regarded as vulnerable because of low allelic richness (Table 4), the small number and size of populations, and the subspecies decline (Lee & Walsh-McGehee, 2000). Various conservation interventions should include control or eradication of invasive predators (especially rats and domestic cats), and restoration of breeding habitats.
Phaethon lepturus ascensionis includes all populations of the southern Atlantic Ocean (i.e., Abrolhos, Fernando de Noronha, Ascension Island, São Tomé, and Príncipe). This subspecies is particularly at risk on the Brazilian islands where populations are small, threatened by introduced predators and exhibiting low genetic diversity and distinctiveness between the tiny yet discrete populations at Abrolhos and Fernando de Noronha (Nunes et al., 2017). Nunes et al. (2017) estimated the effective population size of Fernando de Noronha to be approximately 100 birds. A recent genetic bottleneck was detected in the Ascension population using the heterozygosity excess method (Appendix S1, Table S5) but not using the M-ratio. Such method-related differences have important implications for assessing the timing of population declines with such bottleneck tests (Peery et al., 2012). In principle, the population decline occurred more recently based upon heterozygosity excess estimates. The total population size of this subspecies may be <2,000 Príncipe (50-100 pairs) and São Tomé (unknown number of pairs; Bollen, Matilde, & Barros, 2018). These populations are threatened by introduced predators (cats, rats, and tegus Salvatore merianae; Nunes et al., 2017). In Brazil, white-tailed tropicbirds are listed as "Threatened" on the Brazilian Red List . Because of small population sizes and the low number of breeding colonies in the subspecies P. l. ascensionis, we consider this conservation unit to be threatened.
Due to the distinct genetic nature of P. l. europae, the Europa Island population should be regarded as a single-island endemic unit. This subspecies was described by Le Corre and Jouventin (1999) on the basis of its morphometric measurements and plumage color.
Microsatellites and mtDNA analyses confirm that P. l. europae is genetically isolated from all other populations, and the mtDNA phylogenetic tree also supports the hypothesis that this population may have been founded from a common ancestor shared with the P. l. ascensionis group (Figure 6). Genetic diversity of P. l. europae was low (Table 4) and reduction in population size was detected using the heterozygosity excess tests under the IAM model (Appendix S1, Table S5). Europa was estimated to contain >1,000 pairs of birds in 1974 (Barré & Servan, 1988), 500-1,000 pairs in 1997 (Le Corre & Jouventin, 1997), and probably < 500 pairs in 2017 (Amy, TAAF, MLC personal communication). Therefore, the population is declining rapidly and conservation measures are needed urgently to reverse this trend. The main threat is predation of eggs and chicks by introduced black rats that are abundant on Europa Island (Russell, Ringler, Trombini, & Le Corre, 2011). Since 1995, breeding success of both red-tailed and white-tailed tropicbirds on Europa is consistently <10% and in some years the latter experience entire breeding failure (Le Corre & Jouventin, 1997;Ringler, Russell, & Le Corre, 2015).
With such low productivity, annual recruitment of new adults is likely to be close to zero and this undoubtedly will result in local population extinction in the near future. Restriction to a single-island, rapid population decline and present-day threats from introduced mammals place P. l. europae in the critically endangered category and stress the urgent need for conservation measures.

| Species complex of P. l. lepturus, P. l. dorotheae, and P. l. fulvus
The separation of P. l. lepturus from P. l. dorotheae and P. l. fulvus is more questionable than the Atlantic and Mozambique Channel subspecies. Population genetic structure analysis showed that P. l. fulvus was inseparable from P. l. lepturus although this population contains an unusual color morph that might be the result of phenotypic plasticity (Price, 2006). Moreover, when we consider mtDNA, P. l. fulvus possesses a private haplotype (Table 4 and Figure 3), suggesting an original maternal lineage. However, the phylogenetic tree did not clearly isolate the population as unique phylogroup ( Figure 6). Our results do not support the validity of P. l. fulvus subspecies. The bottleneck test showed no variation in population size (Appendix S1, Table S5) but analysis of a larger sample size is therefore strongly recommended. The population contains an estimated 6,000-12,000 pairs (BirdLife International, 2018) and is not threatened, although introduced predators and historic deforestation may have impacted the species as it has on all endemic avifauna on Christmas Island (James & McAllan, 2014). Even with the validity of P. l. fulvus subspecies not support by our findings, the Christmas Island population with a differentiated genetic pool may at the very least be considered as a distinct conservation unit and it should merit further study.
Data from morphometric comparisons, population genetic structuring, and phylogenetic analyses indicate that there is no distinction between P. l. lepturus and P. l. dorotheae, represented by the Tahiti population, suggesting that the latter may not be a valid subspecies or the Tahiti population may be considered as the same subspecies as birds in east Indian Ocean populations. Although some of these populations may be threatened (especially by invasive predators), the subspecies P. l. lepturus appears not to be at risk currently. They include hundreds of populations, and tens of thousands of breeding pairs. Furthermore, different populations seemed to be genetically connected, implying a functional metapopulation of birds ( Figures 5 and 6). This means that local population declines may be compensated for by immigrations of birds from other populations.
Dispersion in a metapopulation helps to maintain overall viability by rescuing local populations from potential extinction (Akçakaya, Mills, & Doncaster, 2007), also providing hybrid vigor and preventing inbreeding depression (Saccheri & Brakefield, 2002). The mean allelic richness per locus of the studied populations of P. l. lepturus subspecies and the Tahiti population ranged from 5.5 to 6.8, which was higher than the allelic richness of all other conservation units of the species (Table 4) and indicates a potential for adaptive evolution of these colonies. None of these populations appeared to have experienced genetic bottlenecks, except for the Mauritius population (Appendix S1, Table S5), for which a recent founder effect is