Similar to Chen et al. (2016) and Peng et al. (2020), the present study area is located in inshore waters surrounding the DRE in the northern apex of BBG, northern South China Sea ( Figure 1 ). The present datasets were collected by a joint effort of two independent teams from the Institute of Deep-sea Science and Engineering, Chinese Academy of Science (IDSSE), and Shantou University (STU). IDSSE’s surveys followed predetermined zigzag survey routes to ensure consistent coverage across the study area. STU’s surveys were conducted by searching along the coast back and forth at different distances from the coastline. Given the small range of the DRE humpback dolphins, the whole study area could be well covered within 1 day regardless of the field methods adopted. For both teams, regular surveys were carried out year-round using a small boat powered by a 60-HP engine. Dolphins were searched by naked eyes; once encountered, photos of both sides of their dorsal fins and upper bodies were photographed with digital SLR cameras (Canon EOS 7D, Olympus E-M1 Mark II) equipped with either 100–400- or 100–200-mm-zoom lenses, irrespective of their distinctiveness, age, or behavior. Each encounter was followed for a minimum of 10 min or until all group members were covered by the photographer. Geographic locations were recorded using GARMIN 78s.

Photos were processed and filtered using the program DISCOVERY (Gailey and Karczmarski, 2012) according to their qualities (1 ≤ Q ≤ 100) concerning the focus, the contrast from the background, the angle of the animal to the camera, and the proportion of the body region used for identification (Karczmarski et al., 2005). Each dorsal image was scored from 0 (not marked) to 5 (highly marked) according to its distinctiveness (Friday et al., 2000) including the fin shape, notches on the trailing edge of the dorsal fin, pigmentation, and permanent injuries on the upper body region (Karczmarski et al., 2005; Lin et al., 2018). Moreover, the DRE humpback dolphins were categorized into three age groups, the calves, juveniles, and adults ( Table 1 ). Calves were excluded from the following analyses regardless of their distinctiveness as they generally have different demographic dynamics compared to the juveniles/adults.

Mark-recapture models assume no individual heterogeneity in their demographic perspectives; the violation of this assumption may lead to biased estimates of parameters (Sandercock, 2006). Unfortunately, there is no goodness-of-fit test available for the robust design model used in the present study, so the model fit could not be evaluated or adjusted. Instead, over-dispersion of sighting data was assessed by calculating the variance inflation factors ( c ^ ) of a saturated Cormack–Jolly–Seber (CJS) model (Cormack, 1964; Jolly, 1965; Seber, 1965) using the median- c ^ method (Anderson et al., 1994). As a rule of thumb, c ^ < 3 indicates adequate fit to the data and the impact of data over-dispersion is negligible. Test2.CT and Test3.SR embedded in the program U-CARE (Choquet et al., 2009) were also used to test for “trap effect” and “transient impact”, respectively.

We applied the Huggins robust design model framework (Pollock, 1982) embedded in the program MARK (White and Burnham, 1999) to infer the demographic parameters of a subset of the present dataset. The robust design model gives estimates of five parameters, including φ_t (the probability that individuals survive from primary occasion t to t+1), γ′ and γ″ (the probability of being off the study area/unobservable during primary occasion t given that the individual was unobservable or observable at time t-1, respectively, and survives to time t), and p _ts and c _ts (the probability that the individual is captured or recaptured in secondary occasion s of the primary occasion t, respectively).

The robust design consists of two levels of sampling occasions. The first or primary sampling occasions are separated by a relatively long period, and therefore the population can be geographically (through migration) and demographically (through birth/death) open. The primary occasions consist of multiple secondary sampling occasions during which the population is considered closed. The other major assumptions made by robust design include: 1) marked individuals that do not lose their marks and whose marks are not overlooked; 2) no individual heterogeneity in capture/recapture within each sampling session; 3) no individual heterogeneity in survival between primary occasions; and 4) the fate of individuals, including the probabilities of being captured and surviving to the next occasion, being independent.

In the present study, the sighting histories were pooled into secondary sampling occasions every 2 months from 2015 to 2019 (Silva et al., 2009; Chan and Karczmarski, 2017), which should be short enough to consider such a small population closed (Zeng et al., 2020), but sufficient to reach a high capture rate. The bimonthly sighting data with <40 records were discarded as the low sighting rate would result in imprecise estimates of parameters, with the rest bimonthly data collapsed into the year as primary occasions ( Table 2 ). Effort was also made to reduce the bias due to other assumption violations. For assumption 1, we used only the images with high quality (Q ≥ 70) and moderate-to-high distinctiveness (Q ≥ 3) in the subsequent analysis (Friday et al., 2000; Karczmarski et al., 2005), to reduce the bias associated with mismatch. For assumptions 2 and 3, we reduced the difference in survival by excluding neonate/calf data and also divided the other individuals into different age/sex groups which may present different demographic dynamics. Assumption 4 has to be violated for dolphins as they live in social groups, which may lead to overdispersion of the sighting data.

The marked ID-ratio ( θ ^ ) was calculated as the proportion of moderately-to-highly marked dorsal fins (D ≥ 3) over the total captured dorsal fins in each sampling year (excluding the calves), as we assumed no difference between the capture probabilities of marked and unmarked individuals. The total population size (including non-marked, low and highly-marked individuals, N ^ T ) was then reconstructed by correcting the marked population size (N) by the marked ID ratio ( θ ^ ) following: N ^ T = N ^ / θ ^ . The variance of N ^ T was calculated following Urian et al. (2015) as:

The upper and lower bounds of N ^ T were estimated with N ^ T lower = N ^ T / C and N ^ T upper = N ^ T × C , where:

(Burnham et al., 1987)

Site fidelity of the DRE humpback dolphins was measured by calculating the lagged identification rate (LIR), which represents the probability of individuals identified within the study area to be identified again t time units later. The LIR is expected to be constant if the population is closed. A reducing LIR with time lag indicates that the probability of recapturing individuals in the study area drops due to either emigration or mortality (Whitehead, 2001). If the LIR drops but remains at a level above zero or bounces back after a certain time point, it indicates that some individuals reside in or re-immigrate to the study area. Models assuming closed population, emigration or mortality, and emigration + mortality were subsequently fitted to the observed LIR using SOCPROG (Whitehead, 2009), and the 95% CI of model parameters was calculated using 1,000 bootstrap replications. The QAIC criterion was used to rank and identify the optimal models (see Results for details).

During the 5 years of study (2015-2019), a total of 27,915 images were collected through 112 days of field survey. Of these photos, 198 individuals were successfully identified, including 147 highly marked individuals (D ≥3). The cumulative discovery curve sharply increased at the early stage (April 2015–March 2016) of fieldwork but slowed down in the 2nd year and reached a plateau in early 2017 ( Figure 2A ). Seventeen new IDs captured after April 2017 were exclusively young juveniles or calves entering the marked population with emerging identifiable marks. Individual dolphins were sighted an average of 13.9 times (range from 1 to 48, Figure 2B ). Individuals with only one sighting record (n = 22) comprised 12.2% of the dataset, while 123 individuals (83.7%) were sighted more than five times.

U-CARE showed some evidence of the “transient impact” (Test3.SR, χ² = 26.90, p < 0.01), which was primarily contributed by the sex-undetermined adults (χ² = 23.21, p < 0.01) rather than the juveniles (χ² = 5.72, p = 0.06) or female adults (χ² = 2.72, p = 0.44). Test2.CT provided no evidence of “trap effect” (χ² < 0.001, p = 1.00). The signed square root of the χ² statistic (z) was estimated as zero, suggesting that the difference between the capture probabilities between newly captured and recaptured individuals was particularly small. The median c ^ was estimated as 1.68 (SE = 0.25), indicating a limited impact of data over-dispersion. Given the fact that CJS and robust design models are built on different mathematics frameworks, and that the γ of the robust design model should account for the “transient impact,” we did not adjust the robust design modeling results with the c ^ value.

The annual mark ratios ranged from 0.921 to 0.958 during the present study period, with a mean value of 0.934 (SE = 0.005) for the DRE humpback dolphin ( Table 3 ).

The full model received the least statistical support (model #27 in Table 4 ). No evidence was detected for trap effect (model #26 vs. #27, ΔQAIC_C = 40.5) or group-specific capture probability (model #24 vs. #27 , ΔQAIC_C = 68.5). There was little improvement in model fit when the survey effort was incorporated as covariate (model #25 vs. #24, ΔQAIC_C =1.1). Therefore, only time effect on the capture probabilities (p = c(t)) was used for all the subsequent analyses.

During the second round of modeling, the model considering both time and group effects received no statistical support (model #21 in Table 4 ). The model considering a constant survival (φ(.)) best fit the data (model #2). The model in which a constant survival was attributed to female adults and a time-dependent survival was attributed to sex-undetermined adults and juveniles (φ(J=UA(t))/FA(.)), model #4) received a very close statistical support to the one of model #2 (ΔQAIC_C = 0.29), followed by the ones assuming age-specific and constant survivals (φ(J(.)/FA=UA(.)) in model #6, ΔQAIC_C = 1.35 compared to model #2). Based on these three model structures on survival, a series of models were further built to test for seven movement hypotheses. As shown in Table 4 , there was no specific movement hypothesis that best described our data. However, the model fit was generally improved when the movement parameter number was reduced by adding constraints to γ (e.g., even flow), especially when either γ′ or γ″ was set to specific values (e.g., no emigration and no movement).

The averaged estimates of capture probabilities ranged from 0.36 to 0.90, with the mean value as 0.60. Eleven out of the sixteen secondary occasions presented capture probabilities > 0.5 ( Figure 3B ). Survival presented a similar temporal pattern for the juveniles (0.92, SE = 0.05; 0.86, SE = 0.06; 0.88, SE = 0.03; 0.89, SE = 0.03 from 2015 to 2018, respectively) and sex-undetermined adults (0.92, SE = 0.05; 0.86, SE = 0.06; 0.88, SE = 0.03; 0.89, SE = 0.03 from 2015 to 2018, respectively), with a drop in 2016 followed by a slight increase in 2017 and 2018 ( Figure 3A ). A constant survival was recorded for the female adults (0.89, SE = 0.03).

γ was fixed in five out of the eight optimal models (model #2~6). The values averaged over the rest optimal models suggested a small γ″ (γ″₂₀₁₅ = 0.03, SE = 0.04; γ″₂₀₁₆ = 0.08, SE = 0.07; γ″₂₀₁₇ = 0.02, SE = 0.04) and a large γ′ (γ′₂₀₁₆ = 0.88, SE = 0.20; γ″₂₀₁₇ = 0.97, SE = 0.04, Figure 3C ), while the last γ′ and γ″ were not estimable.

After being corrected by the mark ratio, the optimal robust design models suggested an annual population size of 156 (95% CI: 133–184), 133 (95% CI: 120–148), 114 (95% CI: 107–122), 105 (95% CI: 102–108), and 102 (95% CI: 98-107) from 2015 to 2019, respectively ( Figure 4 ). The population change parameter (N_t/N_t+1) varied from 0.853 to 0.973, with a mean estimate of 0.901 (SE = 0.06).

The QAIC criteria suggested the “Emigration + re-immigration” as the optimal model to fit the observed LIRs of the DRE humpback dolphins ( Supplementary Table 1 ). The best-fit model suggested a mean value of 2556.8 days (95% CI = 85.2–8519.8) and 263.8 days (95% CI = 12.7–10993.8) for the rest time in and outside the study area, respectively.

When different age and sex groups were examined separately, the hypothetical model considering mortality and emigration received the best support in fitting the juvenile LIRs ( Supplementary Table 1 and Figure 5 ), with the mean residence estimated as 2,899.3 days (95% CI = 1549.1–9530.0); the QAIC provided comparable support for most of the hypothetical models in fitting the sex-undetermined adults’ LIRs (ΔQAIC_C < 2, Supplementary Table 1 and Figure 5 ), and the mean residency varied from 332.8 days (95% CI = 0.4–8095.1) under the “Emigration + re-immigration” hypothesis to 2,2306.7 (95% CI = 5401.9–8.4 × 10¹³) under the “Emigration/mortality” model. For the female adults, the “Closed” population model and “Emigration + re-immigration + mortality” model received comparable statistical support (ΔQAIC_C = 0.68), with the latter suggesting a rapid turnover of dolphins using the study area (mean time in = 33.2 days, 95% CI = 29.6–632.9; mean time out = 3.2 days, 95% CI = 0.1–116.2).

Two major indexes are used to assess the robustness of the photo-ID dataset, including the population coverage and the capture probability. First, sufficiently high coverage of the population is particularly important for the population size estimation; otherwise, the mark-recapture model will underestimate the population size since part of the population members do not have a chance to be captured. Second, low capture probabilities, which are generally associated with insufficient sampling effort, may reduce the accuracy and precision in the population size and other parameter estimates (Tyne et al., 2016). Lin et al. (2018) simulated a sighting matrix data with a series combination of different population sizes and capture probabilities and showed that capture probabilities <0.5 will lead to >10% bias in the N estimate and capture probabilities <0.3 will lead to over 5% uncertainty (relative standard error) in the parameter estimate for small populations (~100 in population size).

Since the cumulative discovery curves of both of the two previous studies (Chen et al., 2016; Peng et al., 2020) and the present study reached the plateaus, the DRE humpback dolphins were expected to be covered by all the three independent studies. Of these, Chen et al.’s (2016) study reported a highly comparable catalogue size (N = 151) to the present study (N = 147 highly marked individuals), while Peng et al.’s (2020) catalogue was notably larger (N = 230, >50% larger compared to the former two studies). It is noteworthy that these three studies covered the same area in a consecutive time frame from 2011 to 2019. Therefore, the difference in catalogue sizes was unlikely caused by the natural demographic process (e.g., birth or immigration) for two major reasons: 1) long-lived marine mammals like the humpback dolphin cannot experience over 50% expansion through birth within such a short time; and 2) the DRE humpback dolphins were generally considered geographically closed (see Residency/Fidelity), and the transience, if any, should be encompassed by the present study based on a much larger dataset with longer sampling duration (60 months) compared to Peng et al.’s (2020) research (31 months).

Alternatively, artificial error, i.e., the false reidentification of individuals (re-sighted individuals have been added as new members), seems to be the most plausible explanation. The humpback dolphins in Chinese waters are well-known for their skin discoloration process (Jefferson and Leatherwood, 1997), which provides an excess of recognizable features for individual identification. Nevertheless, given that these features may change fast (Guo et al., 2020), it raises the risk of false reidentification, especially when the sighting interval is too long due to an insufficient survey effort (Lin et al., 2018; Guo et al., 2020). False reidentification may lead to a low capture probability and an overestimated population size. Unfortunately, neither Chen et al. (2016) nor Peng et al. (2020) reported their capture probability estimates. Considering that over half of their sighting data was composed of individuals ≤2 sightings (17% of the present dataset for comparison), a bias associated with low capture probabilities and/or insufficient survey effort is probable.

Ns derived from different mathematical frameworks differ in their biological meanings. Specifically, the N of robust design or POPAN represents the number of individuals using the study area during the primary occasion (annual population size in the present study) or throughout the entire sampling period (2011–2014 in Chen et al., 2016), respectively. If the population is demographically constant (birth rate was equal to mortality rate), the N of POPAN should be [number of annual births × (number of sampling years -1)] larger than the N of robust design. If the population was declining through increasing mortalities, the difference in N estimates will further include the increasing number of deaths. Therefore, the Ns of different studies cannot be compared directly (Sandercock, 2006).

As expected, the present study proposed slightly smaller, yet comparable, N_m estimates (annual marked population size: 95–144) in comparison to Chen et al.’s (2016) N_m (total marked population size: 159). However, the difference in N_t estimates (102–156 in the present study vs. 261 in Chen et al. (2016) is too large to be justified by the difference of their biological meanings. Instead, this difference should be primarily due to the mark ratio used for N_m correction. For instance, the mark ratio used by Chen et al. (0.61, 2016) was substantially smaller than the range repeatedly calibrated by most of the other conspecific studies, including the humpback dolphins from the Xiamen Bay (0.94, Zeng et al., 2020), Hong Kong waters (0.86, Chan and Karczmarski, 2017), the east Taiwan Strait (0.86-0.94, Wang et al., 2012), and the present study (0.92–0.96). The cause of the low mark ratio in Chen et al.’s study is unknown. When readjusted with the average mark ratio of the present study, it would result in a total population size (N_t) of 170 (SE = 13, 95% CI = 146–199). If taking this readjusted N_t as the closest, albeit overestimated, estimate of the initial population size during Chen et al.’s (2016) study period in 2011, it would suggest 8.2% of individual loss from 2011 to 2015, which might be further accelerated (34.6% of individual loss) from 2015 to 2019.

The adults’ survival of DRE humpback dolphins was strikingly low. One may argue that the survival was underestimated due to models underestimating the emigration rate, implying that the models failed to distinguish the effect of mortality and emigration. However, independent photo-ID studies reported no individual exchange between the DRE and the neighboring populations since the early 2000s (Chen et al., 2016; Tang et al., 2021), which received support from the genetic differentiation recently confirmed by using microsatellite loci and Hi-C data (Zhang et al., 2020; Zhao et al., 2021). Therefore, long-distance migration, if any, should not affect the demographic process of this population. In other words, the loss of individuals of the DRE humpback dolphins should primarily represent mortality.

The survival of the DRE humpback dolphins was substantially lower than the range reported for conspecific populations inhabiting the highly urbanized Xiamen Bay (0.976, Zeng et al., 2020), west Taiwan coast (0.985, Wang et al., 2012), and Hong Kong waters (0.980, Chan and Karczmarski, 2017). Consistent with this low survival, an abnormally high carcass encounter rate was recorded during our field study. For instance, three fresh carcasses were spotted in a single month in 2015. Although the stranding or mortality information was not systematically collected, there is no doubt that most of the deaths went unobserved. To our knowledge, such low survival was rarely reported for delphinid populations, e.g., spinner dolphin (Stenella longirostris) in Hawaii (0.97), bottlenose dolphins (Tursiops truncatus) in the Azores archipelago (0.97), and Sarasota (0.92-0.99) (Silva et al., 2009; Tyne et al., 2014). One exception is the bottlenose dolphin in Barataria Bay, whose survival dropped from 0.95–0.96 to 0.87 after the Deepwater Horizon oil spill in 2010 (Lane et al., 2015). Unfortunately, the cause(s) of death remained undetermined for the carcasses recovered in the DRE region (Zhu et al., 2019), preventing us from identifying the major threat(s) to the population. However, these potential threats may be considered severe and intense since they caused acute mortality as suggested by the fluctuating survivals.

Survivals differed between the age/sex groups of the DRE humpback dolphins. In general, the juveniles and the male adults present lower survival compared to the female adults (Fruet et al., 2015; Arso Civil et al., 2019), which was generally explained by the size or sex selective pressures. For the DRE humpback dolphins, however, the survival of the juveniles and sex-undetermined adults fluctuated around the value of female adults, suggesting that the dominant pressure not only varied on different age groups but also presented a rapid temporal change.

Of all the candidate threats, that associated with dolphin-watching activities which occur primarily within the administrative boundary of the Qinzhou City is of particular concern (Wu et al., 2020). On the other hand, dolphin ranging patterns may vary between sex and age groups. For instance, female bottlenose dolphins in Shark Bay presented strong locational philopatry (Tsai and Mann, 2013), as the female calving success may benefit from staying in the familiar habitats and associates (Rendell et al., 2019). Female-biased philopatry might also exist in the DRE humpback dolphins (see Residency/Fidelity), which would account well for the relatively constant survivals observed for the female adults in the present study. In comparison, the point estimates of the juvenile and sex-undetermined adult survivals, although with overlapping confidence intervals, seem to imply an improving survival, which was consistent with the reducing population decline rate of the DRE humpback dolphins.

The robust design model incorporates the effect of temporary emigration and therefore is expected to provide a more accurate estimate of survival. However, the residency of female adults, which presented the lowest level of survival, was best described by no movement (“closed population”). Although the robust design modeling detected certain evidence for the temporary emigration, the emigration probability was estimated not higher than 0.08, i.e., only a few individuals. Consequently, it should not provide sufficient supply for the re-immigrants. Thus, the emigration of the present study is likely to be artifact, and the reducing identification/sighting probability of individuals is most likely determined by mortality instead of animal movement.

To date, two studies compared the photo-ID catalogues of the six humpback dolphin populations in China, which are usually separated by over 150 km, and found no match (Wang et al., 2016; Tang et al., 2021). The lack of permanent or temporal migration was also reported in smaller geographic scales. For instance, Guo et al. (2020) found no individual overlap between the two flanking zones of the Pearl River Delta region (ca. 50 km in linear distance). An absence of individual exchange was also reported between DRE humpback dolphins and the neighboring population in its eastern waters (ca. 60 km in linear distance, Chen et al., 2016). These findings, albeit based on relatively short survey periods, support the lack of long-distance travel or high levels of site fidelity and residency of the humpback dolphins in Chinese waters. Therefore, the movement outside the study area as revealed by LIRs is most likely to be explained by the temporary use of the edge of the population range, rather than the migration between populations.

The present study presented a robust photo-ID dataset to assess the movement and demographics of the DRE humpback dolphins. By reviewing and comparing with the previous demographic parameters, it was suggested that the DRE humpback dolphins have experienced a rapid decline in abundance in the mid-2010s, with the latest estimate of 102 (95% CI: 98–107) in 2019. The population remains unsustainable with the present survival rates, especially considering the small size of population. Therefore, it is urgent to identify the major threat(s) and design conservation/management programs accordingly in the near future.

Eco-tourism has been proposed to be one of the major threats to the DRE humpback dolphins (Chen et al., 2016; Wu et al., 2020). Dolphin–boat interaction has been suggested to alter the dolphin behavior budget (Constantine et al., 2004). If the behavior harassment is intense, it may further cause the displacement of animals (Shannon et al., 2017). An ongoing study suggests that the core habitat of the DRE humpback dolphin had shifted from the waters off Sanniang Cove (where it is also the dolphin watching site) during the 2000s to the shallow waters off the Liangjiang Bay in the early 2010s (Lin et al. in prep.). Although the dolphin–boat interaction is generally not considered a lethal factor, the risk associated with living in an un-preferred area may increase. Therefore, the indirect impact of dolphin–boat interaction, especially from a demographic perspective, should be included in conservation plans and dolphin-watching activities should be regulated in the future.