a) BBS dataset
The BBS monitoring program, led by Endemic Species Research Institute in Taiwan, has been conducted since 2009. The aim of the BBS is to monitor the long-term population dynamic of breeding birds in Taiwan. The BBS dataset included 457 BBS sites located across the Taiwan island from 2009 to 2017 (Figure 1), ranging from 0 m a.s.l. to 3900 m a.s.l. Each BBS site included 6 to 10 points located within an area of 2×2km, and each point was spaced at least 200 m apart.
The BBS surveys were conducted by point counts from local sunrise to 4 hours after local sunrise in good weather conditions (i.e., no rain during the survey). The surveyor counted and recorded the number of all the birds heard or seen for six minutes at each point in three distance bands (0–25, 25–100, and >100 m). Birds heard or seen
twice in each year except of the year 2009, which was visited three times in a year. These two visits of a given site should be at least two weeks apart. In order to match the bird breeding season at different altitudes of Taiwan, low-elevation sites (<1000 m a.s.l.) were surveyed once in March and once in May; mid-elevation sites (1000–2500 m a.s.l.) were surveyed once in April and once in June; and high-elevation sites (>2500 m a.s.l.) were surveyed once in May and once in June. Each visit of a BBS site included a total duration of between 36 to 60 minutes (6–10 points) and a total survey area of between 0.1884 to 0.3140 km2 (based on the 100 m radius circles).
Among the 142 BBS sites originally established since 2009, only 27 BBS sites (19%) were continuously surveyed until 2017 (Table S1). The Endemic Species Research Institute recommended that a maximum of four surveyors could participate in each visit, in order to control the effect of number of surveyors on the survey. From 2009 to 2017, only 0.42% of the 4949 visits had five or more surveyors. The average of observed species richness reported from each point was 7.16 species (Figures S1). The average of observed species richness reported from each visit was 15.78 species (Figures S2). In the rank abundance distribution plot of all BBS survey data from 2009 to 2017, I applied the Null, Preemption, Log-normal, and Zipf models evaluated by Bayesian Information Criterion (BIC). The Log-normal distribution (BIC = 43335) has the best fit among all the models for rank abundance distribution from a total of 4949 visits (Figure S3).
Figure 1Location of the 457 BBS sites surveyed on Taiwan island from 2009 to 2017
b) eBird dataset
I downloaded eBird data recorded from 1967 to 2018 through the eBird database (https://ebird.org/data/download/ebd). eBird database has recorded more than four hundred thousand checklists in Taiwan (as of July 2020) (https://ebird.org/taiwan/home).
Four primary survey protocols have been defined in eBird – stationary, traveling, historical, and incidental. The definitions of the four primary survey protocols are as follows: (1) stationary survey protocol follows in a single fixed location with no more than 30 m away from the starting point of the checklist, and the surveyor is required to know the exact start time and duration. According to the eBird’s survey protocol recommendation, duration under three hours makes the better information of the checklist (i.e., shorter checklist gives scientists more accurate information about the exact location and time of birds occurrence); (2) traveling survey protocol follows a distance with more than 30 m away from the starting point of the checklist, and the surveyor is required to know the exact start time and duration. In addition, the specific distance of traveling is required to submit or the surveyor needs to estimate the distance traveled to the best of their ability. The eBird’s survey protocol recommends keeping traveling distance under eight km in order to make a better quality of checklists; (3) historical survey protocol only requires the surveyor to know the date of birding. In other words, the exact time of day, duration, and distance traveled are not required to submit. In some cases, historical checklists may consist of historical bird watching events. For example, data from the Taiwan Bird Record of Chinese Wild Bird Federation (CWBF), had recorded 102,716 checklists from 1972 to 2017 (Lin et al., 2020). However, some locations and duration reported were not accurate from the CWBF dataset; (4) incidental survey protocol refers to those checklists which bird watching is not the primary purpose (e.g., attention might be focused on driving, gardening or doing indoor activities). Incidental checklists lack
important survey information (e.g., the exact start time, duration, and distance traveled) and are less useful for scientific purposes. The eBird database also notifies “complete checklists,” which surveyors report all bird species they were able to detect and identify (does not exclude species or report only highlights). On the other hand, an “incomplete checklist” happens when surveyor intentionally omits any wild bird species that was present, detected, and identified (exclude introduced species, invasive species, and heard or seen-only species). Still, it is feasible to omit any captive species.
In the rank abundance distribution plot of all eBird data recorded in Taiwan from 1967 to 2018, I applied the Null, Preemption, Log-normal, and Zipf models evaluated by BIC. The Log-normal distribution (BIC = 907987) has the best fit among all models for rank abundance distribution from a total of 313,050 eBird checklists (Figure S4). In addition, the three most common sampling protocols each made up nearly one-third of the total dataset: stationary (31.21%), historical (31.34%), and traveling (30.35%), incidental (7.07%) (Figure S5). Checklists with a duration of ≥6 minutes made up 93.48% of the dataset (Figure S5 and Figure S6).
Methods 1. Bird data
a) BBS dataset
I obtained BBS dataset recorded from 2009 to 2017 through the Endemic Species Research Institute, Taiwan (https://sites.google.com/a/birds-tesri.twbbs.org/bbs-taiwan/bbs-zi-liao-shen-qing). I selected data which were recorded from March to July. I excluded BBS sites that contained less than 6 points and BBS data that were recorded farther than 100 m from each point. I only included bird species that regularly breed in Taiwan during the breeding season. A total of 135 diurnal resident and summer visitor bird species from BBS dataset were included in this study (Table S2). Thus, non-breeding bird species (i.e., wintering, transient migrant, pelagic seabird, vagrant, and introduced species) were all excluded throughout the study (Table S2). The migratory statuses of bird species followed the 2020 Checklists of Birds of Taiwan, Chinese Wild Bird Federation.
To make our results comparable to the eBird database, I only selected BBS sites which included at least six completed and approved eBird checklists within a 2x2 km square buffer based on centroid point from each BBS site with ArcGIS 10.6. More than half of the BBS sites (55%) included less than six completed and approved eBird checklists (Figure S7). The main principle for establishing BBS sites is based on the criteria to include national parks, important bird and biodiversity areas (IBA), and wildlife refuges, which represents the complete breeding bird community and environment in a particular area. The BBS sites established along the coast are intended to include more types of habitats (habitat heterogeneity). Thus, to exclude the main habitats of the most wintering, transient migrant, and pelagic seabird species, I removed BBS sites which were intersected with coastline. A total of 204 BBS sites were retained after selection (Figure
2). Among the 204 remaining BBS sites (n = 2238 visits), 165 sites were located in low-elevation (<1000 m a.s.l.); 29 sites in mid-low-elevation (1000–2500 m a.s.l.), and 10 sites in high-elevation (>2500 m a.s.l.) (Table S1).
Figure 2 Distribution of selected 204 BBS sites (orange-colored) across Taiwan island from the original of 457 BBS sites from 2009 to 2017
b) eBird dataset
I included eBird dataset recorded from March to July, 2008 to 2018. As described above, I only included diurnal birds that breed in Taiwan. A total of 144 bird species from eBird dataset were included in the study (Table S2).
I selected the completed and approved checklists which were intersected within a 2x2 km square buffer based on centroid point from each BBS site with ArcGIS 10.6. A total of 2591 locations were reported across Taiwan’s main island (Figure 3). If any location where eBird checklists uploaded was intersected from two or more BBS sites at the same time, I treated eBird checklists separately belonging to each BBS site; though, this rarely occurred.
To avoid duplicate checklists in the eBird and BBS datasets, I excluded eBird checklists with location names that had similar patterns to BBS sites, such as “BBS-A35-19”. For survey protocol selection, I selected checklists from the three most common survey protocols, as follows: stationary, traveling, historical (including data uploaded from the Taiwan Bird Record of Chinese Wild Bird Federation). I only included checklists that were at least 6 minutes in duration for the comparison to the BBS dataset (Figure S6). Based on the two primary high intensity periods of bird activity during a day (Robbins, 1981), I restricted eBird checklist start times to after 4 AM and end times to before 7 PM (Figure S8). The number of surveyors in each eBird checklist was mostly under four persons (Figure S9), which matches the BBS survey protocol of including under four surveyors in each visit.
To minimize misleading results of species richness estimation in subsequent analyses, I removed the whole checklist if any bird species was reported as “X” (no
Incidence-based species richness estimation, which only requires to submit presence-absence data, I transformed any species reporting more than one individual to “1”. I removed species independently with the individual count which obtained “NA” (no data available) from the report. To represent the presence of a species, I transformed any species reporting as “X” to “1”, without removing the whole checklist. Also, I removed the duplicated checklists, which were usually shared by individuals of same birding group, based on the sampling event identifier. Eventually, a total of 14596 checklists that fell within BBS sites were collected for further analyses.
Figure 3 Distribution of eBird checklists reported locations across Taiwan. A total of 2591 locations were reported from 2008 to 2018.
2. Statistical analysis
a) Observed species richness comparison
The different BBS survey methods employed in 2009, caused a different time of duration in each visit than other years (i.e., 6-minute point count surveys were conducted from 2010–2017, while 9-minute point count surveys were conducted in 2009). I therefore removed all visits from the BBS dataset from 2009.
To make the results comparable, I compiled species records and duration of survey points of a given BBS site in a visit. After compiling records of a visit into a checklist in each site separately, a total of 2238 checklists were collected from each visit across the 204 BBS sites in Taiwan. To be comparable with BBS’s survey duration, I only included eBird checklists with a duration of between 36 to 60 minutes, with a total of 2164 eBird checklists retained. I performed a two-tailed Wilcoxon rank-sum test on both datasets to test the difference of observed species richness.
b) Species richness estimation methods
For the selected 14596 eBird checklists that fell within BBS sites, species richness estimation was based on each separate checklist (checklist-based). Three non-parametric approaches of species richness estimation methods were applied to the eBird dataset: (1) abundance-based estimator, Chao1 (Chao, 1984; Colwell & Coddington, 1994; Chao &
Chiu, 2014); (2) Incidence-based Coverage Estimator (ICE) (Chao & Chiu, 2014):
recommended by Chao and Chiu (2014), I set up 10 individuals as a cut-off point to define infrequent or frequent species group; (3) and first-order Jackknife, an estimator based on the number of singleton species (Burnham & Overton, 1978; Colwell & Coddington, 1994). Chao1 estimation was performed using the “iNEXT” package (Hsieh et al., 2016);
ICE and first-order Jackknife estimation methods were performed with the “vegan”
package (Oksanen et al., 2016) in the R platform.
c) Evaluating the performance of species richness estimation methods
To quantify the performance of the species richness estimation methods from the eBird dataset, I calculated the bias value based on estimated species richness from each eBird checklist against the compiled observed species richness from 2009–2017 in each BBS site separately (i.e., the asymptote of total species richness from accumulated annual surveys was assumed to be known as the total species richness in each BBS site, likely to represent the local bird community) (Walther & Morand, 1998; Walther & Martin, 2001;
Walther & Moore, 2005; Tingley et al., 2020). In other words, each eBird checklist produced one result value of bias (unless the eBird location was intersected with more than two BBS sites, then I treated the eBird checklists separately belonging to the shared BBS sites). The bias value was calculated by the following formula:
Bias = [𝑬𝒊𝒋#𝑨𝒊] [𝑨𝒊]
with j = eBird checklists in the i th BBS site (i.e., j th sample in each BBS site);
with i =1 to 204 (refers to the i th BBS site). Eij is the estimated species richness in each eBird checklist; Ai is the compiled observed species richness of the i th BBS site from 2009 to 2017. The bias calculation was performed in Microsoft Excel 2019. Finally, I used one-tailed Wilcoxon rank-sum test to examine the least biased species richness estimator among the three estimation methods by comparing each pair of estimators. The selected least biased species richness estimator was applied to the species richness estimation in order to access the two datasets comparison in the following questions.
d) Determining the effect of duration on bias after species richness estimation (1) Evaluating the effect of duration on observed species richness
Before taking the next step to examine the effect of duration on bias, I tested the effect of duration on observed species richness across all included 14596 eBird checklists.
I fitted four non-linear functions independently by using the least squares method (James et al., 2013). The four non-linear functions are used to estimate the asymptote of species richness as duration increase (Magurran & McGill, 2011), and formulas are depicted as follows:
where, y is the observed species richness, as the dependent variable, and x is the duration, as the independent variable; a, b, c denote the parameters to be estimated by the least squares method. This parameter estimation was calculated with the “stats” package (Team & Worldwide, 2002) in the R platform.
To compare the goodness-of-fit of the four different non-linear models, I compared the fitted curve with the BIC (Gideon, 1978). BIC was used instead of Akaike information criterion (AIC), since our objective was to explain the relationship between duration and observed species richness, instead of predicting the value (Shmueli, 2010).
Under the Bayesian probability framework, the probability of selecting the true model increases as the training sample size increases (Friedman et al., 2001; Magurran & McGill, 2011). BIC model selection was performed with the “AICcmodavg” package (Mazerolle
& Mazerolle, 2019) in the R platform. The best selected non-linear function was used to address the relationship between the duration and bias in the following process.
(2) Calculating the reduction of bias after species richness estimation
To make a comparison of the reduction of bias before and after estimating species richness at a standardized duration, for the same reasons as above, I removed all visits from the BBS dataset from 2009. With a total of 14596 eBird checklists, I treated duration in each eBird checklist as an independent variable; bias derived from observed and estimated species richness were treated as a dependent variable separately. Bias was calculated by the following formula: eBird checklist (note that the estimation was based on the least biased estimation method);
Ai is the compiled observed species richness from the i th BBS site recorded from 2010 to 2017.
To test the effect of duration on the bias across all included 14596 eBird checklists, I fitted both independent and dependent variables with the selected non-linear function described above by using the least squares method (James et al., 2013). Parameter estimation was calculated with “stats” package (Team & Worldwide, 2002) in the R platform. Finally, based on the non-linear function at a 60-minutes, the reduction value of bias can be measured with – the bias value after species richness estimation minus the bias value before species richness estimation.
(3) Evaluating improvement on proportion of species richness from eBird against BBS after species richness estimation
To evaluate the improvement of species richness after estimation from eBird dataset against BBS dataset at the duration of 60 minutes, I included BBS sites which only included 10 points (i.e., a total of 60 minutes in each visit was conducted from a BBS site), and removed all visits from 2009. I calculated the average observed species richness from each visit in each BBS site (i.e., the average number of species recorded in each visit of BBS). A total of 92 BBS sites were retained after selection (Figure 4), accompanied with a total of 6611 eBird checklists. I treated duration in each eBird checklist as an independent variable; bias derived from observed and estimated species richness were treated as a dependent variable separately. Bias was calculated by the following formula:
Bias = [𝑶𝒊𝒋#𝑨𝒊] [𝑨𝒊]
with j = eBird checklists in the i th BBS site (i.e., j th sample in each BBS site);
with i =1 to 92 (refers to the i th BBS site). Oij is the observed species richness in each eBird checklist; Ai is the average observed species richness from each visit in the i th BBS site recorded from 2010 to 2017.
Bias = [𝑬𝒊𝒋#𝑨𝒊] [𝑨𝒊]
with j = eBird checklists in the i th BBS site (i.e., j th sample in each BBS site);
with i =1 to 92 (refers to the i th BBS site). Eij is the estimated species richness in each eBird checklist (note that the estimation was based on the least biased estimation method);
Ai is the average observed species richness from each visit in the i th BBS site recorded from 2010 to 2017.
To test the effect of duration on the bias across all included 6611 eBird checklists, I fitted both independent and dependent variables with the selected non-linear function described above by using the least squares method (James et al., 2013). To test the performance of eBird dataset after species richness estimation, based on the non-linear function, 60-minutes was set to standardize the comparison of bias before and after species richness estimation. Finally, the improvement on proportion of species richness from eBird dataset after the estimation can be calculated through the bias formula.
Figure 4 Distribution of selected 92 BBS sites with criteria of 10 points/site from 2010 to 2017 across Taiwan from the original of 457 BBS sites
Results 1. Observed species richness
After restricting duration from both BBS and eBird datasets (36–60 minutes/checklist), the BBS dataset (204 sites) had a statistically higher observed species richness than the 2164 eBird checklists which were recorded within a 2×2 km square buffer based on centroid point from the BBS sites (W = 3826200, effect size = 0.503, p <
0.001) (Figure 5). The median per checklist of observed species richness for BBS (n = 2238) and eBird (n = 2164) datasets were 15 and 9 species, respectively. Inter-quartile range (IQR) for BBS (n = 2238) and eBird (n = 2164) datasets were 9 and 8, respectively (Figure 5).
Figure 5 Observed species richness per checklist recorded in BBS and eBird datasets.
BBS dataset included 2238 visit-based checklists, with a total of 204 sites. eBird dataset included 2164 checklists. Both datasets had durations restricted to the range of 36–60
2. The performance of species richness estimation methods
Chao1 estimator (median bias = -0.693) was overall least biased (W = 12369000, p < 0.05) compared with other two estimators (median bias of ICE = -0.730; median bias of Jackknife = -0.773) against compiled observed species richness from each BBS site (Table 1, Table 2 and Figure 6). ICE estimator was less biased than Jackknife (W = 119220000, p < 0.001) (Table 2). Estimates of species richness by eBird checklists varied by estimation methods, but generally underestimated the true community size (bias < 0) (n = 14596) (Table 1). However, the outcome of estimated species richness varied across estimation methods. Bias derived from the Chao1 estimator varied between -0.987 and 5.602, while bias derived from the Jackknife estimator has a generally smaller range, varied between -1.000 and 1.000 (Table 1).
Table 1 Performance of three species richness estimation methods for the eBird dataset against observed species richness from the BBS dataset, evaluated by the result value of bias summarized by all included checklists (n = 14596). Bias was calculated to make a comparison among estimators.
Mean SD Median IQR Minimum Maximum
Chao1 -0.576 0.393 -0.693 0.440 -0.987 5.602
ICE -0.640 0.286 -0.730 0.351 -0.983 1.222
Jackknife -0.689 0.267 -0.773 0.317 -1.000 1.000
Table 2 One-tailed Wilcoxon rank-sum test between species richness estimation methods
W–value p–value
Chao1 vs. ICE 123690000 < 0.05*
Chao1 vs. Jackknife 123690000 < 0.05*
ICE vs. Jackknife 119220000 < 0.001***
Figure 6 Performance of Chao1, ICE, and Jackknife estimators on species richness estimation methods. Bias was measured by comparing the result of each estimation method against compiled species richness from each BBS site. (A) The difference of Chao1 subtracted from ICE estimator; (B) The difference of Chao1 subtracted from Jackknife estimator; (C) The difference of ICE subtracted from Chao1 estimator; (D) The difference of ICE subtracted from Jackknife estimator; (E) The difference of Jackknife subtracted from Chao1 estimator; (F) The difference of Jackknife subtracted from ICE estimator. Asterisks in plots indicate the significance level between estimation methods by one-tailed Wilcoxon rank-sum test (p < 0.05 = *; p < 0.001 = ***). Note that the result value of bias only presents from -0.05 to 0.05.