All fieldwork was conducted according to permits N° 5343/2016, delivered by “Servicio Agrícola y Ganadero” (SAG). The handling of birds was kept to the shortest possible time, and a hood was systematically used to reduce stress.

The samples used in this study were collected during two field campaigns (austral summer: January 2015 and Austral winter: July 2015). Work was conducted on Motu Nui, a small islet located 1.3 km off the Rano Kao crater in Rapa Nui, Chile (27° 20’S, 109° 45’W). Detailed sample sizes are given in Supplementary Material Table S1 and Supplementary Materials Figure 1 .

Assigning an individual bird to a particular species in the Pterodroma genus is not a trivial achievement (Brooke and Rowe, 1996; Brown et al., 2010; Brown et al., 2011). As a prior for species assignment, we used the taxonomic classification proposed by Onley and Scofield (2016). The criteria for species identification relied mostly on jizz and plumage colouration, principally the colour of primary and secondary feathers and underpart and upperpart coverts (Plaza et al., 2021). Several individuals with unclear phenotype or lacking sufficient photographic material for posterior confirmation were classified as “unknown Pterodroma species”. However, as opposed to what is observed in the Round Island hybrid swarm (Brown et al., 2010), we did not observe any individual with a clear hybrid phenotype. We retained only those individuals that could reliably be assigned to one species (N=81). For each individual, we collected photographic documentation, and identification was confirmed independently based on this material by seabird identification expert Peter Harrison (see in particular Harrison et al., 2021 for the identification criteria used). These identifications are referred to here as a-priori assignments, as they do not rely on the objective measurements performed in this study but rather on overall appearance and expert knowledge.

For all individuals, we measured eight morphological traits ( Table 1 ). Since our primary interest is interspecific discrimination, morphometry was transformed into a set of allometric traits: all measures except for body mass were log-transformed and scaled by log-transformed body mass. Body mass itself was log-transformed only. We verified the homogeneity of variance in the transformed traits between species using Levene’s test, accounting for multiple testing through the Benjamini-Yekutieli (BY) procedure (Benjamini and Yekutieli, 2001). We first performed an analysis of variance (ANOVA) to detect differences in size between the six a-priori species, including molecular sex and its interaction with the species, as additional categories and assessed pairwise species effects using Tukey’s test of honest significant differences. Next, we performed a multivariate analysis of variance (MANOVA) on the log-transformed measurements scaled to log-transformed body mass, including a priori species and sex as predictor variables, to detect allometric differences between species. In all tests, we revised the residual distributions to validate the models.

To get a clearer picture of species differences in multivariate morphospace, we performed a linear discriminant analysis (LDA) on the log-transformed measurements scaled to log-transformed body mass, together with the log-transformed body mass itself, using R v. 4.0.3 (R Core Team, 2020 - this version is used throughout) and R package MASS v. 7.3-53.1 (Venables and Ripley, 2002). We reduced the dimension to the first two axes of variance based on the high redundancy of morphometric variables (see Results) - while the choice of LDA (maximising separation between groups), as opposed to PCA (maximising total variance), was guided by the aim of this study. We computed the standard ellipse (SE) for each species in the plane defined by the first two discriminant functions using the Bayesian approach implemented in the R package SIBER v. 2.1.5 (Jackson et al., 2011). Note that although SIBER was initially designed to analyse isotopic data, our data’s multivariate normality respected all the SIBER model’s assumptions for SE estimation - and it offered a more robust approach to SE uncertainty with small sample sizes than analytical approaches. We used “uninformative” conjugate normal priors (µ=0, σ² = 10³) for the group means and Inverse-Wishart (2x2 scale matrix with 2 degrees of freedom) for the covariance matrix and ran two MCMC chains for 10,000,000 generations, with 1,000,000 generations burn-in and retaining every 100th sample. We computed the distribution of overlaps between SEs for each species pair as indicative of the objective morphometric similarity between species.

We used the ratio of nitrogen and carbon stable isotopes (δ¹⁵N and δ¹³C) as a proxy of trophic niche overlap between species. Analyses were conducted both with feathers and whole blood. Feathers are a slow-growing tissue, renewed only during the annual moult, and are expected to integrate stable isotope intake during their growth period, which mostly encompasses the non-breeding period. On the other hand, blood is a fast-turnover tissue that is expected to reflect stable isotope uptake in the weeks preceding the sampling (Quillfeldt et al., 2005; Cherel et al., 2008). Stable Isotope ratios were analysed as described in Zapata-Hernández et al. (2014): Briefly, samples were combusted in a continuous helium flow elemental analyser, and resulting gases were reduced and dried. Finally, the remaining N₂ and CO₂ gases were separated by gas chromatography and analysed by mass spectrometry (Delta Plus spectrometer, Thermo-Finnigan). We first checked the correlation between feather and blood isotopic ratios for individuals for which both samples were taken using ordinary least-square regression and f-test. We then compared δ¹³C and δ¹⁵N separately amongst species using ANOVA (treating feather and blood values separately). We included species, molecular sex (and their interaction), and collection season as categorical variables for feathers (all blood samples were collected in winter). Note that we could not obtain sufficient blood for isotopic analysis in P. nigripennis, so it was excluded from blood analyses. Further, we also compared the position of each species in the bivariate (C, N) plane using the Bayesian ellipse approach implemented in the SIBER package in R. Algorithm parameters were set as above.

In order to assess species-specific patterns in the spatial structure of nest distribution, we georeferenced all Pterodroma nests in both breeding seasons (austral winter or summer, i.e. June or January survey), and used K-means clustering based on pairwise Euclidean distances to define the optimal number of clusters in the nest distribution using the NbClust v. 3.0 package in R, computing all available indices and choosing based on consensus (Charrad et al., 2014). Finally, differences in species composition between clusters were investigated using the Cluster package v. 3.0.1 in R (Maechler et al., 2022). Then, to assess potential differences in nesting habitat between species, we used a nest classification scheme that considered both structure and materials ( Supplementary Material Table S2 ). Six types of nests were observed: (i) exposed, i.e. laying directly on a flat, usually rock surface (ii) inside a rock crevice, (iii) protected by a rock (e.g. against a ledge), (iv) protected by vegetation, mostly long grass (typically grass of the species Paspalum forsterianum, DuBois et al., 2013), (v) protected by a combination of rock and grass, often against the edge of large rocks in the vegetation, and (vi) entirely underground burrow excavated in the dirt (see Plaza et al., 2021). Relative usage frequencies of each nest type by each species were used to calculate niche breadth using Levins’ standardised index (a value of 1 indicating generalist and 0 for complete specialisation). We also calculated niche overlap between all pairs of species Horn’s niche overlap index (with smaller values representing lower overlap). Both indexes (Levins and Horn) improved species’ niche measures and resource preferences (Krebs, 1998).

In order to establish the in-flight vocal repertoires of each species, we recorded aerial calls during the evening hours as birds flew back to the colony. Calls were recorded using a TASCAM DR-07 MKII portable recorder (with a frequency range of 44-48,000 Hz) and a Sony ECM-672 microphone (frequency range of 50-16,000 Hz). Recordings were discarded when the vocalising individual could not be identified with high certainty or when calls had too much background noise. Calls were digitised with a sampling frequency of 24,000 Hz (16-bit resolution) using Raven Pro (Charif et al., 2010). In the spectrograms, we measured the first ten notes of calls based on two principal characters: 1) fundamental frequency or first harmonic and 2) entire note or bandwidth. Three frequency measures were included for each character: the low, high, and peak frequency (HZ). Temporal features were measured with note length (see Supplementary Materials Table S3 and Figure S2 ). Principal component analysis (PCA) was used to obtain a lower-dimensional representation of these bioacoustic variables using the dudi.pca implementation in R. We used PCA scores for linear Discriminant Analysis (LDA) to maximise group separation.

We sampled blood from the brachial vein of breeding birds that we captured by hand directly in the nests. Blood samples were stored on FTA™ Classic Cards, and extracted using a standard NaCl protocol (Aljanabi and Martinez, 1997). We sequenced a 782-bp fragment of the mitochondrial Cytochrome-b gene (Cyt-b) using primers designed from Pterodroma sequences obtained in GenBank (P. nigripennis U74343.1, P. ultima HQ420379.1, and P. neglecta U74341.1): gadpe_1F: GTT CTA CTG GTT GGC TGC CTA C, gadpe_1R: ATG GCC CCA AAC CTT CGA AAR TC, and internal primer gadpe_C22: ATA CG ATA CCR AGA GGG TTG. PCR conditions consisted of 3 min initial denaturation at 95°C, followed by 35 cycles of 30-sec denaturation at 95 °C, 30-sec annealing at 58°C, and 1 min elongation at 72°C, followed by a final elongation step at 72° for 5 min. PCR product purification and bi-directional sequencing were performed by Macrogen (Seoul, South Korea) using the same PCR primers. We also used 45 primer pairs to amplify a panel of microsatellite loci: 9 designed explicitly for Pterodroma petrels (Brown and Jordan, 2009) and 36 universal primers for birds (Dawson et al., 2010). Each locus was tested for polymorphism and fluorescently-labeled using a forward primer modified with an M13 tail, as described in Schuelke (2000). We followed a touch-down PCR protocol, with an initial denaturation at 95°C for 1 min, followed by 27 “hot” cycles (94°C denaturation for 30 sec, 59°C annealing for 30 sec, and 72°C elongation for 30 sec), and 8 “cool” cycles (94°C denaturation for 30 sec, 53°C annealing for 30 sec, and 72°C elongation for 30 sec), and 7 min final elongation at 72°C. PCR products were checked on a 1% agarose gel with a ladder of 500 pb, and analysed on an ABI PRISM 3500 XL (Applied Biosystems) microcapillary electrophoresis analyser. Genotypes were scored using GeneMaker v1.85 (see Supplementary Material Table S4 for details).

Sequences were processed, aligned, and trimmed using Geneious v. 2021.1 (^© 2005-2021 Biomatters Ltd). We calculated haplotype (Nei and Tajima, 1981) and nucleotide (Nei and Li, 1979) diversity, as well as general statistics, using R packages ape v. 5.5 (Paradis and Schliep, 2019) and pegas v. 1.1 (Paradis et al., 2004; Paradis, 2010). In order to gain a clearer view of phylogenetic relationships in the context of Pterodroma diversity, we downloaded published sequences for 13 other Pterodroma species or subspecies (P. arminjoniana, P. deserta, P. feae, P. hasitata, P. incerta, P. macroptera macroptera, P. macroptera gouldi, P. madeira, P. magentae, P. mollis mollis, P. phaeopygia, P. sandwichensis, and P. solandri), previously published sequences for P. alba, P. heraldica, P. neglecta, P. nigripennis, and P. ultima, as well as an outgroup, Pelecanoides garnotii. We used all available sequences if less than six, and chose six random sequences if more, with a total of 183 sequences, including those in the present study. All Genbank accession numbers and references are provided in Supplementary Material Table S5 . We then reconstructed the cytochrome-b gene tree using a Bayesian framework implemented in BEAST2 v. 2.6.4 (Bouckaert et al., 2014). Cyt-b sequences were partitioned into first and second codon positions on the one hand, and third codon positions on the other, specifying a separate (but identical) HKY substitution model for each partition, and estimating both the transition/transversion rate and the underlying nucleotide frequencies, with a gamma-distributed per-site mutation rate discretised into four classes. Considering we were working at an interspecific level, we set a Yule model as a prior on tree topology. MCMC was run for 10,000,000 generations, discarding the first 20% as burn-in and recording every 1,000th state, ensuring the collection of large effective sample sizes (ESS > 200 for all parameters) to estimate posterior parameter distributions. The reconstruction was repeated three times with identical parameterisation to assess convergence further. Although 15 individuals for which species could not be reliably assessed were removed from all other analyses, they were included here to completely inform the fine-scale tree topology (e.g., by providing intermediate sequence states). However, these sequences were never attributed to a particular species and remained in the tree as unspecified sequences.

We calculated expected and observed heterozygosity, and population-level summary statistics for each locus using the R packages adegenet v.2.1.3, pegas, and hierfstat v. 0.5-9. Pairwise Fst was calculated according to Weir and Cockerham (1984), and significance levels were assessed against 1,000 random permutations of individuals among the pairs and adjusted using the BY procedure. We performed admixture analysis using STRUCTURE 2.3.4 (Pritchard et al., 2000) without prior knowledge of species grouping. We used the admixture version of the model, assuming correlated allele frequencies. The MCMC was run with a burn-in of 10,000 steps, and 1,000,000 sampling steps, with many possible clusters (K) ranging from 1 to 10. Analyses were repeated ten times. We used Structure Harvester (Earl, 2012) and Evanno’s delta-K (Evanno et al., 2005) for model choice and summarised the results across the ten independent runs using CLUMPP v. 1.1.2 (Jakobsson and Rosenberg, 2007). In parallel, we used the discriminant analysis of principal components approach (DAPC, Jombart et al., 2010, which does not rely on the Hardy-Weinberg equilibrium assumption) using the adegenet R package (Jombart, 2008). We retained the first 25 principal components of the eigen-decomposed microsatellite variance-covariance matrix, effectively retaining > 99% of the total variance for the linear discriminant analysis. Linear discriminant functions were computed both (1) between K-means-clustered groups (with many groups ranging from 2 to 6), mirroring the Structure approach, and (2) using a priori species as a classification factor. Finally, we used the snapclust algorithm (Beugin et al., 2018) as a third, independent control, using AIC model choice to identify the optimal number of clusters.

We tested the association between genetic and ecological distances for all quantifiable axes of variation at the interindividual and interspecific levels. We used Procrust tests implemented in the vegan v. 2.5-7 package (Oksanen et al., 2020) to correlate distance matrices. At the individual level, we used (a) genetic distance based on microsatellites, (b) genetic distance based on mtDNA Cyt-b sequences, (c) distance in allometric morphospace, (d) distance in feather-based isotopic niche, and (e) distance in the blood-based isotopic niche. At the species level, we used (a) Fst based on microsatellite data, (b) standard ellipse overlap in principal LDA components in allometric morphospace, (c) standard ellipse overlap in feather isotopic niche, (d) standard ellipse overlap in blood isotopic niche, and (e) standard ellipse overlap in principal LDA components in acoustic space. We did not include mtDNA sequences at the species level since the haploid, non-recombining nature of these markers makes them less suitable than nuclear markers when comparing groups (Teske et al., 2018), and we considered microsatellite markers, in that case, adequately captured genetic distance. Acoustic recordings were not associated with particular individuals and could only be used at the interspecific level. Finally, for each pair of species, we scored the number of dimensions along which they could be distinguished, considering (a) phenology, (b) nest choice, (c) isotopic niche, (d) morphology, (d) acoustics, and (e) Cyt-b haplotype. This count was compared to pairwise microsatellite-based Fst between the two species. The correlation between ranks in the differentiation order and Fst values was assessed using Spearman’s rank test.

Based on Levene’s test, log-transformed morphological measurements did not significantly deviate from the homoscedastic hypothesis (all p-values > 0.05). Log-transformed body mass differed significantly between species, but not between sexes, even taking into consideration the interaction of species and sex (ANOVA Dfs, F-values and p-values respectively: [5, 28.5, < 0.001], [1, 0.104, 0.748] and [5, 0.651, 0.662] - see Table 1 ). The differences between P. nigripennis and the other species accounted for most of this observation, as P. nigripennis was, by far, the smallest species (Tukey’s post-hoc test, all pairwise comparisons involving P. nigripennis p-values < 0.001). A significant ordering in body mass could be observed with P. neglecta > P. heraldica > P. nigripennis (all p-values < 0.05). Other pairwise tests were not significant (see Supplementary Information Table S6 ).

The MANOVA analysis showed that species determined a priori was a significant determinant of log-transformed morphological measurements (p-value < 0.001), indicating allometric differences between species (see Supplementary Information Table S6 ). Maximising variance between a priori-determined species in multivariate space through linear discriminant analysis yielded similar observations ( Figure 1 ). The first two discriminant functions accounted for 88.5% of the total between-group variance (75.9% and 12.6%, respectively). The first function’s main loading was wing length (94.2%), and discriminated mainly P. nigripennis from the five other species. The second function’s main loadings were wing length and bill length (39.1% and 54.4%, respectively), allowing most to separate P. ultima from the five other species. The third function covaried with bill length, tarsus width, and wing length (56.0%, 20.4%, and 17.2%, respectively) but accounted for only 7.2% of the total variance. Standard ellipses fitted to species projections in LDA space were consistent with the other results ( Supplementary Figure S3 ): P. nigripennis had no overlap with any other species, and the proportion of overlap was highest between P. alba, P. atrata, and P. heraldica. P. ultima showed low overlap with P. alba, P. atrata, and P. heraldica, but higher with P. neglecta.

Overall, we found that P. nigripennis was separated from all other species in the morphospace, primarily due to its smaller size and relatively longer wings. Based on bill length, P. ultima was distinguishable from all other species, except P. neglecta. The other species could be differentiated pairwise but presented a considerable overlap in their trait size distribution when considered collectively. P. heraldica, P. alba, and P. atrata were nearly indistinguishable based on morphometrics alone.

Stable isotope ratios (δ¹⁵N and δ¹³C) in blood and feathers were loosely correlated when all species were taken together (but significantly with f-test for linear regression at the α = 0.05 level: R^{2 =} 0.10, p = 0.022 for δ¹⁵N and R^{2 =} 0.26, p < 0.001 for δ¹³C - see Supplementary Figure S4 ). Interspecific differences mainly drove the stronger correlation for δ¹³C; consistently with this, analyses of variance revealed that species was the best explanatory variable for δ¹³C, considered a proxy for spatial distribution, and δ¹⁵N, a proxy for trophic level, both in blood and in feathers. Within species, no correlation between blood and feather values was observed for either element, and sex differences were never significant. Carbon signature segregated clearly P. ultima from all other species when analysed in feathers (P. ultima has a significantly lower δ¹³C); when analysed in blood, P. ultima had a significantly lower δ¹³C than all other species except for P. nigripennis, and a significantly ranking could be observed, with P. ultima < P. neglecta < P. heraldica ( Figure 2 ). Note that in the Southern Pacific, lower δ¹³C values generally correspond to higher latitudes (Magozzi et al., 2017). Nitrogen signature segregated P. alba and P. heraldica from all other assessed species when measured in blood, and P. atrata had a lower signature than all other species except P. ultima. The results were similar when measured in feathers, although contrasts were less clear (P. heraldica had a higher ratio than all species except for P. alba, and the lowest-ratio species, P. nigripennis, had a significantly lower ratio than the two highest-ratio species P. alba and P. heraldica). The standard ellipse approach ( Supplementary Material Figure S5 ) showed that isotopic niche width (as approximated by standard ellipse area) was usually higher in feathers than in blood - which may reflect wider-ranging and less specific foraging in the non-breeding period. As expected, the overlap was lowest between P. ultima and all other species.

Pterodroma species on Motu Nui are partly segregated along the year: P. nigripennis was only found in summer (January), while P. alba, P. heraldica, and P. ultima were only observed in winter (June). P. atrata and P. neglecta were encountered throughout the year. Nests appeared grouped into 3 clusters of multispecific composition ( Supplementary Material Figure S1 ), both in summer and in winter: in other words, Pterodroma petrels tended to cluster together with other Pterodroma petrels regardless of their species, as is usually observed for this group throughout their Pacific range (P. Harrison, pers. obs.). Each species preferentially used a different microhabitat in summer to build the nests (chi-square test, p-value = 0.001). P. nigripennis, in particular, was highly specialised in burrow nests, in keeping with observations throughout its range. Similarly, P. neglecta preferred rock crevices, as observed in several localities throughout its range (P. Harrison, pers. obs.). There was little overlap with other species ( Tables 2 , 3 ). In winter, no significant differences were found in substrate nest preference between species (p-value = 0.63). Nest use was less specific, and species overlapped broadly ( Tables 2 , 3 ). Rock-crevice and rock-protected nests were most frequent, and exposed nests were less used.

Aerial calls were recorded around sunset and tended to fade as night passed. P. ultima was never engaged in vocal activity during our fieldwork, and could not be recorded for this study. We recorded a succession of chattering “ti-ti-ti” calls for P. alba, P. atrata, P. heraldica, and P. nigripennis (see Supplementary Material Figure S2 ), with syllable duration that ranged from 0.03 s in P. atrata to 0.11 s in P. nigripennis ( Supplementary Material Table S3 ). The highest fundamental frequency was observed in P. nigripennis at 2178 kHz ( Supplementary Material Table S3 ). These four species seemed to diverge most on syllable duration and frequency of harmonics. On the other hand, P. neglecta used a “kuk-u-er” or “moaning call” type of call ( Supplementary Material Figure S2 ). This call was characterised by a low number of elements, a long syllable duration (0.72 s), and a maximum fundamental frequency of 1.270 kHz ( Supplementary Material Table S3 ), rising abruptly in pitch at the end of the call. In a linear discriminant analysis based on species assignment, P. neglecta and P. nigripennis were the most clearly separated groups on the first two discriminant functions ( Figure 3 ).

We sequenced a total of 84 Cyt-b sequences, 782 bp-long, of which 67 were assigned to one of the six a priori species in the study, and 17 originated from individuals whose species could not be reliably determined. We found 17 haplotypes: P. nigripennis and P. ultima had only private haplotypes; P. atrata, P. alba, P. heraldica, and P. neglecta were mostly mixed ( Supplementary Material Table S7 ). Across all sequences, haplotype and nucleotide diversity were high (H = 0.90 ± 0.015 and π = 0.019 ± 0.0093, respectively). Both parameters were of the same order of magnitude for all species, except for P. nigripennis (only one haplotype, see Supplementary Material Table S7 ). We found 77 variable sites: 29 segregates for P. nigripennis and 12 for P. ultima. None segregated in the other species, although P. atrata and P. neglecta had one private mutation (see Supplementary Material Table S7 ). Our cytochrome-b phylogenetic tree reconstruction, which included previously published data, supported the monophyly of most known Pterodroma species ( Figure 4 ): P. arminjoniana, P. hasitata, P. incerta, P. macroptera macroptera, P. macroptera gouldi, P. madeira, P. magentae, P. mollis mollis, and P. solandri were all inferred to be monophyletic; and P. deserta, P. feae, P. phaeopygia, and P. sandwichensis mostly were, despite minor incongruencies. Consistently with other analyses, P. ultima and P. nigripennis were both inferred to be monophyletic too, and sequences produced in other studies clustered unambiguously with our sequences for both species. Additionally, both P. nigripennis ( Figure 4 , taxon 5) and P. ultima ( Figure 4 , taxon 4) were separated from other species in our study by several intermediate taxa. P. atrata ( Figure 4 , taxon 3) tended to group separately from P. heraldica, P. neglecta, and P. alba; however, that taxon also included 2 P. heraldica samples besides 9 of our 11 P. atrata samples; and 2 P. atrata samples grouped elsewhere in the tree. Finally, all P. alba, P. heraldica, and P. neglecta samples (both in our study and published sequences) grouped in several paraphyletic taxa ( Figure 4 , taxa 1 and 2), forming one larger monophyletic clade that also included all P. atrata samples. This paraphyletic group was composed of two sub-groups: the first contained all P. alba samples, together with mostly P. heraldica and P. neglecta, and the second was composed of P. heraldica and P. neglecta. No other systematic pattern emerged in the relationships of these species, as established through Cyt-b sequences. Interestingly, we also underline the fact that the Trindade petrel P. arminjoniana, which is central in the Round Island hybrid swarm where it is known to hybridise with P. neglecta and P. heraldica, is the closest taxon to P. neglecta, P. heraldica, P. alba and P. atrata in our mtDNA tree.

Nine polymorphic loci were retained out of the 45 tested primer pairs (locus characteristics are given in Supplementary Material Table S8 ), which amplified successfully for 96 individuals, among which 79 could reliably be assigned to one a priori species based on phenotype. All loci were found in Hardy-Weinberg equilibrium when species were taken separately, but only one locus (TG03002) when considering all individuals as one population (see Supplementary Material Table S8 ). As expected amongst species, most pairwise Fst was generally high and significant ( Supplementary Material Figure S6 ): yet remarkably, the three comparisons involving P. alba on the one hand, and P. atrata, P. heraldica, and P. neglecta on the other were not significantly different from random, and Fst values between P. atrata, P. heraldica, and P. neglecta, although significant, were relatively low. Clustering approaches yielded similar results ( Figure 5 ), whether they were based on Bayesian or Maximum-likelihood reconstruction of Hardy-Weinberg-equilibrium population-specific allele frequencies (Structure and snapclust), or on linear discriminant analysis (DAPC). In the latter case, the first linear discriminant function was driven by four main alleles at four different loci, making up 50% of allele contributions. No single locus was over-represented (maximum locus contribution: ~25%). The second linear discriminant was more locus-driven, with one microsatellite (TG04061) making up ~78% of total contributions with its five alleles, which were the five largest contributing alleles. All three approaches (Structure, DAPC and snapclust) converged to an optimal clustering in 5 groups. Overall, the Structure and snapclust approaches gave a more nuanced view, and the DAPC approach gave more clear-cut assignments. Both approaches concurred to make P. nigripennis and P. ultima well-defined groups, P. atrata a rather distinct group despite some inconsistencies, P. heraldica and P. alba a similarly composed entity, and P. neglecta a mixed group with its own component as well as a significant heraldica-alba component and some P. atrata admixture components. In some cases, individuals seemed more similar to another species than their phenotypic assignment: this is the case of one P. ultima (identified as P. atrata), three P. atrata (identified as P. heraldica-alba), and two P. heraldica (identified as P. atrata for one method, and P. neglecta for the other).

We found a strong association between genetic distance and most other distance measures at the interindividual level. Microsatellite distance correlated strongly with allometric distances (p-value = 0.001), with feather-isotopes (p-value = 0.003) and with blood-isotopes (p-value = 0.001). Cytochrome-b distances correlated strongly with allometric distances (p-value = 0.003), and more weakly or not at all with blood and feather isotopes (p-values = 0.012 and 0.455, respectively). At the interspecific level, Fst values correlated with feather-based and blood-based isotopic niche overlap (p-value = 0.034 and 0.025, respectively) but not with morphology (p-value = 0.089) or acoustic overlap (p-value = 0.45). In other words, segregation in the foraging niche appeared to be the strongest and most consistent correlate to reproductive isolation in these six groups. Finally, the number of dimensions in which two species segregated was strongly associated with the Fst of that species pair (Spearman’s test: Rho = 0.76, p-value = 0.001), with, on the one hand, P. heraldica and P. alba virtually indistinguishable on both axes, and, on the other hand, P. nigripennis easily separated from all other species ( Figure 6 ).

In keeping with previous descriptions, we found the black-winged petrel (P. nigripennis) to be smaller, with a more slender bill than all other Pterodroma species in this study. Isotopic signals in feathers revealed a relatively lower trophic level, and the species is well identified by its call, breeding phenology, and choice of the burrow that matches earlier descriptions (West and Nilsson, 1994; Hutton and Priddel, 2002). Genetic evidence reveals a deeply monophyletic group within the Pterodroma genus without detectable gene flow with other taxa. Thus, our findings confirm previous phylogenies based on morphology (Imber, 1985) and genetic markers (Nunn and Stanley, 1998; Penhallurick and Wink, 2004).

Although of a less clearly divergent phenotype compared to P. nigripennis, Murphy’s petrel (P. ultima) is well identified by both ecological and genetic traits. It is unambiguously monophyletic in the mtDNA phylogeny, and well distinguished using microsatellite markers. While identified in morphospace, P. ultima is clearly set apart by its particular δ¹³C signature in blood and feathers. This is best explained by year-round spatial segregation at sea from all other Rapa Nui Pterodroma petrels, in keeping with direct evidence obtained on Pitcairn islands: P. ultima forages in oligotrophic oceanic provinces along the subtropical front and the Humboldt current during the breeding period, and as far as the North Pacific Gyre in the inter-breeding period (Clay et al., 2017). Foraging segregation is well known to be a primary barrier to gene flow in seabirds (Friesen et al., 2007a; Friesen, 2015), so the convergence of genetic isolation and marked foraging segregation between P. ultima and its congeners on Motu Nui is consistent with its status as a different species.

The relationships between the Phoenix petrel (P. alba), the Henderson petrel (P. atrata), the Herald petrel (P. heraldica), and the Kermadec petrel (P. neglecta) are all but clear-cut in our study system. These species overlap in nearly all measured parameters, whether genetic or phenotypic, and sorting is often inconsistent depending on the character considered - although we underline the fact that the relatively low sample size for P. alba (5 individuals) may limit our ability to detect finer interspecific differences. P. alba, P. atrata, and P. heraldica are virtually indistinguishable in morphospace, while P. neglecta is only marginally distinguishable as being slightly larger, with comparatively slightly shorter wings. Based on mtDNA Cyt-b sequences, P. alba, P. heraldica, and P. neglecta are indistinguishable, while P. atrata is fairly distinct. Based on microsatellite markers, P. alba and P. heraldica are indistinguishable, P. neglecta is highly variable and often indistinguishable, and P. atrata is generally distinct, although some individuals are similar to P. heraldica. P. neglecta tends to breed more in summer than in winter, and P. atrata more in winter than in summer, but not systematically. P. neglecta also tend to specialise in rock crevices for nesting more than the three other species. Moreover, while these four groups may usually be differentiated visually by an experienced observer, the traits used mainly rely on plumage colouration: P. atrata is unmistakably dark, while P. heraldica and P. neglecta can reliably be told apart only based on the colour of the primary, secondary and tail feather shaft (black in P. heraldica, light-coloured in P. neglecta). P. alba has overall more contrasted colouration patterns.

For all species except P. nigripennis, microsatellite and mtDNA analyses revealed some inconsistencies between molecular and phenotypic assignment. In our sample, one P. ultima individual was consistently linked with P. atrata in microsatellite analyses but with the other P. ultima specimens in mtDNA analysis. Three P. atrata individuals appeared closer to P. heraldica in microsatellite analyses: one unambiguously grouped with other P. atrata sequences in the mtDNA tree, but the two others with the mixed P. heraldica - P. neglecta - P. alba complex. Conversely, two P. heraldica individuals grouped with P. atrata in mtDNA analyses: one also did in microsatellite analyses but not the other. In all these cases, photographic material confirmed the phenotypic identification without ambiguity, and in no case did we observe an individual with a clearly hybrid phenotype, as is the case in Round Island P. arminjoniana, P. neglecta, and P. heraldica (Brown et al., 2010; Brown et al., 2011). Trivial errors (e.g., in tube labelling), though a possible explanation in some cases, are a priori excluded since (1) species identification was the focus of this study, and utmost care was taken in that respect, and (2) the mismatch between mitochondrial and nuclear analysis from the same DNA extract in half of the cases cannot be accounted for in this way. On the other hand, a moderate level of gene flow between species, and a fair amount of phenotypic variation may well account for this observation. The recurrence of atypical colouration convergence between usually distinct groups has recently been conclusively documented in Eudyptes penguins (Frugone et al., 2018), supporting this hypothesis. It should be noted that overall, the Pterodroma species of Motu Nui are far more similar in genetic, isotopic, or morphometric traits than in general appearance or acoustic traits. Thus, a significant amount of gene flow between groups may occur without apparent effects on phenotypic segregation: we suggest that we may obtain similar results in other locations of the Pacific where these closely-related species co-exist.

Among the four less-clearly sorted species (i.e., P. alba, P. atrata, P. heraldica, and P. neglecta), P. alba and P. heraldica are the most similar by far: they were not segregated by any of the quantitative parameters surveyed in the present study ( Figure 6 ). Both present a similar morphology, δ¹³C, and δ¹⁵N signatures, breeding and nesting habits, and were undistinguishable based on our genetic data. This finding converges with previous results based on morphology and juvenile plumage (Murphy and Pennoyer, 1952), mtDNA (Brooke and Rowe, 1996), and microsatellite data (Booth Jones et al., 2017) in other locations. These two groups remain distinguishable by an experienced observer, based on jizz and plumage colouration in particular (Harrison et al., 2021), and some genetic or epigenetic differentiation is undoubtedly involved in maintaining this contrast. However, we consider that there is now enough independent evidence to question the full species status of P. alba: and we suggest that the species may safely be redefined as a subspecies of P. heraldica, named P. heraldica alba.

Although bioacoustic evidence revealed a clear sorting, neither morphology, foraging behaviour, nor genetic evidence supported a binary separation between the two species. One possibility is that acoustic separation is insufficient to act as a leak-proof prezygotic barrier in sympatry, as is usually the case in birds (Bretagnolle, 1996), and the two species are undergoing active hybridisation after secondary contact on Motu Nui, blurring all genetic separation. We did not detect any hybrid phenotype; however, several birds with unclear phenotypes (classified here as Pterodroma sp. or “unknown”) fell within the haplogroups associated with these two species. These are not represented in our bioacoustic dataset since only birds clearly identifiable in flight were included: thus, sampling bias is likely to under-represent hybrids. These two species are also the most similar in morphology, as the primary, secondary and tail feather shaft is the only minimal diagnostic feature – although the typical P. neglecta’s heavier, “bull-necked” jizz distinguishes it from other groups to the experienced observer, and the species includes a characteristic white morph never found in P. heraldica (Harrison et al., 2021). Finally, microsatellite data show a continuous gradient of differentiation, rather than a clear-cut separation, between P. heraldica and P. neglecta. Overall, our results converge with observations made in another recently-colonised secondary contact zone, Round Island, in the tropical Indian Ocean, where admixture between P. heraldica and P. neglecta petrels has been documented by Brown et al. (2010); Brown et al. (2011) and more recently by Booth Jones et al. (2017) – although in Motu Nui, unlike in Round Island, we did not identify any hybrid individual. Thus, the two species do not appear to have strong pre- or post-zygotic isolation, and P. heraldica may be more appropriately described as a subspecies of P. neglecta named P. neglecta heraldica (in that case, we should also consider the current P. alba as P. neglecta alba).

The contrast between P. atrata and the P. heraldica - P. alba - P. neglecta group is easier to define. In the earlier classifications (e.g., Murphy and Pennoyer, 1952), P. heraldica was considered a subspecies of P. arminjoniana (as P. a. heraldica) and further divided into two morphs (“light” and “dark”, Murphy and Pennoyer, 1952), each restricted to specific islands of the Pacific. Imber (1985) first proposed to consider P. heraldica as a separate species from P. arminjoniana. Finally, Brooke and Rowe (1996) proposed considering the dark morph of the new P. heraldica as another species, P. atrata, based on mtDNA divergence between the morphs on Henderson Island, and the then-assumed endemism of the dark morph to that location. Generally, our results support the view that P. heraldica and P. atrata are somewhat distinct (Fst is moderately high, and the morphs are generally distinguishable both on the mtDNA tree and in the microsatellite, structure analyses). However, either ongoing gene flow, or incomplete lineage sorting as a result of very recent reproductive isolation, clearly blur that signal. While we do not see reasons to question the species status of P. atrata based on our data, we find grounds for a more in-depth revision of the phylogenetic relationships between the two species across their range.

This study documents a particular case of sympatric breeding in closely-related Pterodroma species. Seabird species sharing the same reproductive habitat classically segregate in at least one dimension of their niche, such as their at-sea foraging behaviour (Cody, 1973; Young et al., 2010b; Navarro et al., 2015; Ravache et al., 2020). On the other hand, our results document a significant overlap between individuals assigned to different species by the current taxonomy. In particular, overlap in trophic position is striking, and in clear contrast with all other examples to our knowledge, as trophic segregation between sympatric seabird species appeared to be the rule in tropical and subtropical islands (Cherel et al., 2008; Young et al., 2010a; Ravache et al., 2020), and was so far only reported lacking in the subantarctic islands of Patagonia (where the overlap is explained by the high marine productivity and the overall abundance of prey, Forero et al., 2004). This surprisingly promiscuous nesting, however, is not restricted to Motu Nui but seems typical of Pterodroma petrels nesting sites in the South Pacific (Harrison et al., 2021). The extreme local conditions may be the primary driver of nest choice, as proposed by Schreiber (2002): the South Pacific islands share a particularly harsh tropical environment, with intense radiation all year round, high heat, high precipitation, and intense winds. In this context, nest selection may respond more to climatic constraints than to selection for coloniality with conspecifics, and high sympatry, even at very fine local scales, may be expected to arise as a result. We suggest that this immediate daily contact may challenge the weak prezygotic isolation in these closely-related species, and may lead to increased gene flow across species boundaries. Indeed, it may be that high philopatry (and usual allopatry between species) make accurate recognition of conspecifics rather superfluous in most cases for birds that breed in such a fragmented habitat: the recent recolonisation of Motu Nui may pose an unusual challenge to Pterodroma petrels. Such gene flow at secondary contact zones has been evidenced in Laridae (Sternkopf et al., 2010; Sonsthagen et al., 2016), and patterns compatible with hybridisation have been documented in P. arminjoniana, P. heraldica, and P. neglecta in their contact zone in the Indian Ocean (Brown et al., 2010; Brown et al., 2011).

The well-documented local extinction or extirpation events in prehistoric and historical times at the margins of Polynesia make it plausible that our observations are the result of recent secondary contact, and support the idea of Motu Nui as a “melting pot” where gene flow is enhanced between species that are re-extending their ranges. However, these patterns may instead be the result of interrupted ecological speciation (Nosil et al., 2009): indeed, the high overlap of Pterodroma petrels of the alba, atrata, heraldica, and neglecta types may suggest that speciation processes have not led to strong reproductive isolation at any point at all, and that our observations on Motu Nui are just the expected result of the encounter of groups that have never really developed into separate species at all (sensu Mayr, 1996). The positive association between phenotypic and niche divergence, and genetic differentiation, is classically described in that context (Langerhans et al., 2007; Nosil and Sandoval, 2008), and our results unambiguously fit with that model.

Overall, the multidisciplinary approach that we implemented in the present study proves to be a necessary framework to understand speciation processes in closely related species with long generation times, where genetic patterns may be hardly understandable on their own, and we advocate its use in a wider-scale re-evaluation of the taxonomy of the Pterodroma genus. However, the surprisingly high similarity between these currently recognised species would justify the use of higher-accuracy techniques in morphometry, bioacoustics, and in particular in genetics (e.g., genome-wide markers), which may help resolve the finer patterns of admixture in heraldica-alba petrels, or between heraldica and neglecta individuals. Finally, including other populations of Pterodroma species and long-term monitoring spanning several generations is likely to bring key information about the ongoing reversal of speciation – a process of high relevance in the context of massive biodiversity loss.

We underline that whether these patterns result from incomplete speciation or secondary hybridisation, the seabirds of Motu Nui and similar Pacific islets currently constitute an assemblage whose properties are poorly described by the number of species alone. Biodiversity is often summarised through beta-diversity (after the Whittaker framework, 1960) as the ratio of the regional and local species diversity: but we show here that species co-occurrence may involve the appearance of novel diversity. In other words, the true diversity in such areas may not be fully captured by the sum of co-occurring species alone.

Our results also acutely raise the question of the role of species in conservation. While most conservation plans are built on lower taxonomic levels, the treatment of communities with active hybridisation is unclear: these may either represent a threat to the genetic integrity of the species involved or rather a healthy expected aspect of evolution and adaptation. Gene flow across species or subspecies borders seems to be a typical aspect of the Pterodroma genus (Brown et al., 2010; Brown et al., 2011; Booth Jones et al., 2017), and may in itself promote local adaptation in these seabirds’ changing environment. In any case, only long-term monitoring of such communities would allow us to detect possibly positive or negative effects of this cross-species gene flow.

Thus, while questioning species boundaries on Motu Nui, we do not suggest revising the conservation status of the Pterodroma groups involved nor consider ongoing hybridisation as a negative phenomenon. Globally, the P. alba taxonomic unit is considered a species in danger of extinction (BirdLife International, 2021), and should remain so: whether as a subspecies of P. heraldica or as a full species, the threat remains. On the other hand, Pterodroma petrels are collectively defined as a priority action target in the Agreement for the Conservation of Albatrosses and Petrels (ACAP), based on the scarcity of available information about this genus (ACAP, 2017). We strongly recommend that newly arisen assemblages and hybrid groups such as the one on Motu Nui should be closely monitored to evaluate the need for a specific hybrid conservation unit (Piett et al., 2015).