The Role of Urbanness, Vegetation Structure, and Scale in Shaping Puerto Rico’s Acoustically Active Mangrove Fauna Communities

Introduction

Urban bird communities have repeatedly been shown to change in composition with increasing urbanization (Emlen, 1974; Tilghman, 1987; Clergeau et al., 1998, 2001; Chace and Walsh, 2006; Ortega-Álvarez and MacGregor-Fors, 2009). Generally, urban avifauna are most represented by generalist species capable of foraging on a broad range of food types, with the more specialist species largely absent from the most urban sites (Chace and Walsh, 2006). Most of these studies also acknowledge the importance of various metrics of the surrounding land cover matrix, such as forest patch size or habitat diversity, as well as vegetation structure. Further, there is some indication that scale plays an important role in these communities and the relationship between alpha and beta diversity in urban systems (Norris, 2012; Meffert and Dziock, 2013). The above patterns in urban avifauna cover a broad range of habitats and latitudes, although comparatively little is understood about the fauna communities in urban mangroves.

A number of studies have examined urban mangrove avifaunal communities, but there is no agreement on what constitutes “urban” and few conclude on the overall influence of urbanization on bird community structure (Ward, 1968; Lim and Sodhi, 2004; Mestre et al., 2007; Acevedo and Aide, 2008; Mohd-Azlan and Lawes, 2011; Pawar, 2011; Mancini et al., 2018). Several studies acknowledge a relatively high presence of mangrove specialists in anthropogenic habitats of Malaysia, suggesting these species are particularly adept at colonizing urban landscapes (Ward, 1968; Noske, 1995; Lim and Sodhi, 2004). Yet, in Darwin Australia, most of the mangrove specialists have yet to colonize the urban mangroves (Noske, 1996), while mangrove bird richness in general is most dependent upon patch size and the diversity in the surrounding land use matrix (Mohd-Azlan and Lawes, 2011). In São Sebastião, Brazil, differences in mangrove avian communities of a port area are largely attributed to variability in habitats and vegetation structure, with higher species richness associated with greater habitat heterogeneity (Mancini et al., 2018).

In Puerto Rico, studies have shown that urban avifauna communities are distinct from their forested counterparts and that certain traits and distribution statuses (e.g., introduced, endemic, and nocturnal, etc.) are associated with higher or lower occurrences within the island’s urban matrix (Ruiz-Jaén and Aide, 2006; Acevedo and Restrepo, 2008; Suarez-Rubio and Thomlinson, 2009; Herrera-Montes and Aide, 2011). Endemics, for example, seem particularly drawn to the largest forest fragments (Acevedo and Restrepo, 2008; Suarez-Rubio and Thomlinson, 2009), while introduced species, some generalists (Quiscalus niger and Turdus plumbeus), and many anurans (Eleutherodactylus coqui and Leptodactylus albilabris) seem capable of inhabiting urban areas (Ruiz-Jaén and Aide, 2006; Suarez-Rubio and Thomlinson, 2009). For mangroves, Acevedo and Aide (2008) found Twenty-six species of birds in one forest, with forest structural metrics (e.g., canopy cover and tree size, etc.) and the presence of standing water being the most significant variables associated with these communities, although no conclusions are drawn on the influence of urbanization.

Although urbanization is not a primary driver of global mangrove deforestation, it is the most important cause of loss in some cities and is one of the few fragmentation drivers that is permanent (Bryan-Brown et al., 2020; Goldberg et al., 2020). Thus, urban development not only replaces critical habitat for some species, but remaining urban mangrove forests are subject to long-term disturbance forces that create novel systems and faunal communities that may be adapting to these anthropogenic landscapes (Luniak, 2004; McDonnell and Hahs, 2015). Understanding these communities and the underlying forces driving observed changes is then critical to sustainable mangrove management in the Anthropocene. While much of our understanding of avifaunal and anuran communities comes from visual point-count surveys, the present study used passive acoustic monitoring (PAM) to quantify species presence. Acoustic monitoring has repeatedly been identified as being similar to or more effective than visual point counts in a variety of ecosystems (Haselmayer and Quinn, 2000; Zwart et al., 2014; Alquezar and Machado, 2015), and especially in mangroves (Celis-Murillo et al., 2012). Here, we use acoustic surveys and subsequent species identification models to detect avian and anuran presence in twenty forested mangrove sites spanning an urbanization gradient in Puerto Rico. We then aim to answer how mangrove avian and anuran communities are related to various metrics of surrounding urbanization and land cover at different scales, focusing on specific feeding guilds, distribution groups, and mangrove specialists. Understanding these patterns and the potential underlying causes will be critical to managing urban mangrove forests toward optimal function and ecosystem services provisioning.

Materials and Methods

Site Description

The urban mangroves of Puerto Rico, and the sites corresponding to the present study, are described in detail in Branoff and Martinuzzi (2020). Forest composition and structure does vary among sites, but in general the forests are composed primarily of Laguncularia racemosa, followed by Rhizophora mangle and Avicennia germinans. Additional non-halophytic flora is present, and the overall floral richness increases with human population density surrounding each forest. The floral richness of the three mangrove species, however, decreases as urban cover increases. Forest height varies between ten and twenty meters, stem density between five and eight thousand stems per hectare, and basal area between 27 and 43 m² ha^–1. Most of the forest structural metrics are correlated with surrounding vegetation cover or surface water chemistry, with the largest contiguous forests being characterized as taller and with greater variation in height.

Relatively little is known of the mangrove fauna of Puerto Rico in comparison to other ecosystems on the island. The most urban mangroves in the Caño Martĩn Peńa (MPNMAX in this study) were surveyed in 1899 and noted the presence of Megaceryle alcyon, Pelecanus occidentalis, Ardea spp., Rallus sp., Scolopacidae sp., Setophaga spp., Tyrannus dominicensis, and Vireo altiloquus (Evermann et al., 1902). Another more recent survey found the aquatic bird community in the city of San Juan to be composed mostly of Fulica caribea, Bubulcus ibis, Pluvialis squatarola, Pelecanus occidentalis, and Calidris alba, although this study was not limited to mangroves (Fidalgo-De Souza, 2009). Acevedo and Aide (2008) surveyed the bird community of the Levittown mangroves and found twenty-six species, with Coereba flaveola, Quiscalus niger, Parkesia noveboracensis, and Melanerpes portoricensis to be among the most abundant. Descriptions of the mangrove bird community in the islands second largest metropolitan area of Ponce are scarce, but other studies along the southern coast of the island note the use of mangroves for both foraging and nighttime roosting habitat, with the most abundant species being Setophaga petechia and Parkesia noveboracensis (Smith et al., 2008; Rodríguez-Colón, 2012).

Study sites were selected to provide the greatest range in urbanness, at both the island and local scales (Figure 1 and Table 1). To do so, an urban index was developed for all mangrove shoreline across the island, as described in detail by Branoff (2020b). This unitless index constitutes a range of values from 100 representing the most urban site to 1 representing the least urban site. High urban index values are those representing the most urban cover, the least vegetated cover, the highest street density, and the highest population density, and vice-a-versa for the lowest index values. These descriptive variables were pulled from circular buffers surrounding each site at distances of 0.5, 1, and 5 km, and were provided by a variety of spatial land cover and social datasets (Table 2). Three watersheds were first selected as those with mangroves exhibiting the greatest range in urban index values. They are the San Juan Bay Estuary (SJBE) and the Rio de la Plata (Levittown) along the northern coast, and the Rio Inabón to the Rio Loco (Ponce) on the southern coast. Nine locations were selected among these three watersheds, with a minimum and maximum urban site at each. Within the SJBE, seven locations are represented by each of the estuary’s seven waterbodies, with two sites at each location representing the minimum and maximum urban index values for that waterbody. This resulted in 14 sites for the SJBE. Sites are named for their waterbodies and for their local urbanness levels (BAH is the San Juan Bay, MPN, and MPD are the non-dredged and dredged portions, respectively, of the Caño Martín Peña, SAN is the San José lagoon, SUA is the Suárez canal, TOR is the Torrecilla lagoon, and PIN is the Piñones lagoon) and MIN and MAX are the minimum and maximum urbanness sites in each of those waterbodies (BAHMAX, BAHMIN, TORMAX, and TORMIN, etc.). For Levittown (LEV) and Ponce (PON), relatively little mangrove area was available and therefore sites were selected to represent the minimum (MIN), mid (MID), and maximum (MAX) urban index values across the watershed, as opposed to individual waterbodies as was done in the SJBE. This resulted in three sites each in the two locations of Levittown and Ponce. Together, with the seven locations in the SJBE and one each in Levittown and Ponce, a grand total of nine locations and twenty sites are represented by this study.

Figure 1. Passive acoustic monitoring devices were placed in twenty one-hectare mangrove forests spanning an urbanization gradient across the island of Puerto Rico, as represented by the urban index. This index is calculated using the surrounding cover of urban, vegetated, open water, and mangrove classifications, as well as road and population density. Individual locations have a “MAX” and “MIN” urbanness site, referring to the relative urban index values. Urban characteristics for each site are described in Table 1.

Table 1. Site abbreviations used throughout the study and their corresponding locations as well as select urban variables from a 1,000 m sampling radius, as adapted from Branoff (2020b).

Table 2. Covariates used to test hypotheses on the most important drivers of fauna community richness and structure.

Faunal Surveys

Two surveys were conducted in each site, one from March 11th to March 22nd, 2017 and another from November 1st to November 12th, 2017. Hurricanes Irma and María passed near and over Puerto Rico in September of 2017 and changed the forest structure between the two sampling periods, with no available data describing the post-hurricane forests. Therefore, data from both surveys are reported but hypothesis testing is conducted primarily on surveys from the undisturbed pre-hurricane forests. Passive acoustic monitoring devices were fixed to a tree at a height of approximately 1.5 m at each site and programmed to sample for 1 min at 10-min intervals at a sampling rate of 44.1 kHz. Recorders are comprised of an Android smartphone attached to an external Monoprice microphone with flat response between 50 Hz and 20 kHz and a sensitivity of −45 ± 2 dB, all enclosed in a waterproof case. Recorders were programmed with the ARBIMON Touch application and allowed to record until the smartphone battery died, or until they were collected 10 days later. Each site was sampled with an audio recorder device for ∼1 week in each survey period generating ∼1,500 1-min recordings per site/survey.

Acoustic Pattern Matching

All recordings (32,345 1-min recordings) were stored, processed, and analyzed using the RFCx-ARBIMON II platform¹ (Aide et al., 2013). Recordings were visually inspected for interference or damage caused by temporary microphone malfunction or rain and these recordings (6.3% of all recordings, Supplementary Table 1) were removed from the analyses. Species specific acoustic models were constructed by first finding the most representative, loud, and clear examples (templates) of each species’ most consistent vocalization(s). Species identifications were checked with those found on xeno-canto.com and by the authors’ personal records. In the case of Ardea spp. and Psittacidae spp., similar species could not be differentiated by the authors and are referred to by their genus and family, respectively. Call templates were used in a pattern matching algorithm to scan all relevant (diurnal, nocturnal, or both) recordings and identify potential regions of interest (ROIs) with matching patterns. All ROIs with values above the selected correlation threshold (0.1) were presented as potential detections for posterior validations (LeBien et al., 2020). The selection of a low threshold resulted in a high number of false positives though the number of false negatives was negligible. We manually reviewed the pattern matching results. In this step, we annotated the results as either positive or negative, indicating the corresponding species presence or absence. Then, the same templates, along with validated presence and absence recordings were used to train Random Forest Models to recognize the selected acoustic patterns (Breiman, 2001; Aide et al., 2013, Campos-Cerqueira and Aide, 2016). Models were considered sufficiently accurate if they achieved at least 80% accuracy and precision in at least ten of the validated recordings. These models were then used to again scan all recordings for the presence of each vocalization and each species. All recordings marked as present for a given species were again visually scanned by the authors to verify its presence. The presence of each species at each site for each deployment was then cross-checked between the two methods, pattern matching and Random Forest Models. Although most of the cross-validation resulted in agreement, any disagreement between the two was resolved manually by again visually inspecting the inconsistent recordings. Thus, presence for a species at a site for each deployment is defined when at least two of three methods (i.e., patten matching, Random Forest Models, and manual inspection) resulted in presence. For some species, an insufficient number of verified presence recordings were identified for robust Random Forest Models so only verified pattern matching results were used. Detection frequencies were calculated as the percent of surveyed days in which a species was detected.

Statistical Testing

Tests on species richness, defined as the number of species, were performed on groups defined by distribution status, mangrove specialization, and feeding guilds. Avian species were classified into feeding guilds of granivores, omnivores, invertivores, frugivores and nectarivores, and carnivores, piscivores and scavengers following Wilman et al. (2014). Avian endemism was determined from Lepage (2018) and mangrove specialization was determined through (Luther and Greenberg, 2009), which distinguishes between dependent and restricted mangrove specialists. Resident or migrant status was determined through Oberle (2018) and Del Hoyo et al. (1992). Anuran diets were determined from Oliveira et al. (2017) and endemism from Heatwole and MacKenzie (1967).

Variability in species richness was used in hypothesis tests against several covariates describing surrounding land cover, forest composition and structure, and hydrology (Table 2). Surrounding land cover was sampled from a 2 m resolution 2010 dataset of Puerto Rico (Office for Coastal Management, 2017), using circular buffers or radii 0.5, 1, and 5 km. The percentage of each land cover class was calculated as the area of each class within the circle over the area of the circular buffer, and this value was used as a predictor in regression models. Forest composition and structure variables were taken from Branoff and Martinuzzi (2018), and describe the forests before the passing of hurricanes Irma and Maria between the two deployments. These describe summary variables from the 10 plots within each site, and include tree size (diameter at breast height, dbh) and canopy height, as well as species composition descriptors. Finally, site hydrology was described by variables taken from water level models from a 5 years period that were developed with input from water-level measurements recorded at the sites and from regional weather stations (Branoff, 2020b). Hydrology variables include the average water depth, the flooding frequency, and the percentage of time a site is flooded.

All statistical tests were performed in the R language (R Core Team, 2020) and all scripts and raw data can be found in the repository on the Open Science Framework project page at: https://osf.io/uk2fd/. Relationships of detection frequency and faunal richness, with predictors of land cover, forest structure and composition, and hydrology, were tested through simple (one independent variable) linear and ln-linear regression in the lm function. Only the most explanatory predictor variables are reported on and these were chosen as those with the highest R² value among statistically significant models (p < 0.05). Student t-tests were performed through the t.test function in base R to test for differences in species richness between the most urban sites and the least urban sites, and this was done for island-wide urbanness as well as local urbanness. Island-wide urbanness was determined by splitting the sites into the “most urban” and “least urban” by the highest and lowest urban indices, respectively. Local urbanness was determined by the MIN and MAX subscripts at each location, which was also set by the urban index as described in section “Site Description”. For local urbanness, the two MID sites (PONMID and LEVMID) were dropped from the analysis in order to maintain only two categorical levels. Species-specific tests for associations between presence and local MAX and MIN urbanness were performed through Fisher’s Exact Test for Count Data on contingency tables for locations where a species was found only at one of the two sites at a location. This was performed through the fisher.test function in base R. To test for differences in the avian and anuran community structures among the most and least urban mangroves, and to account for differences in call-frequencies, nonmetric multidimensional scaling (NMDS) using the Bray-Curtis dissimilarity index (Bray and Curtis, 1957) was also performed. This was done using the proportion of recordings in which each species was detected at each site, as recommend by Minchin (1987) and as performed through functions of the vegan package (Oksanen et al., 2018). The metaMDS function iterates the NMDS process from random starting points and chooses the solution with the overall lowest stress and thus the most accurate representation of the sites in the ordination space (Clarke and Ainsworth, 1993). To show how site characteristics and species’ detection frequencies correlate with the grouping of the faunal communities, vectors representing site environmental metrics and species were fit to the ordination through the envfit function. Ordination ellipses representing the 95% standard error confidence ellipse of the group centroid locations were drawn for the island-wide least urban and most urban groups as well as for the local MIN and MAX urban groups using the ordiellipse function, also from vegan. To test for differences in communities within the ordination, a permutation multivariate analysis of variance using distance matrices (ADONIS) was performed through the adonis function of vegan. This was also performed on the percentage of recordings at each site in which each species was detected and grouped according to both the island-wide and local relative urbanness characterizations as described above (detection frequency∼urbanness).

Species accumulation curves for each site and each watershed were constructed by counting the number of time-ordered recordings necessary to detect a new species. The specaccum function from the vegan package was also used to construct curves from random permutations of the data and a jackknife extrapolation. Logarithmic curves were then fit to these species-accumulation plots (Ugland et al., 2003), and the number of recordings required to reach a given number of species, as well as the number of species predicted for a given sample area were calculated through the predict function. All graphical visualizations were produced from the ggplot2 package in R (Wickham, 2009).

Results

Five anuran species from two families and thirty-one bird species from nineteen families were detected across all sites during the study (Supplementary Table 1). All five anurans and five bird species are endemic to the Puerto Rico Bank. An additional five bird species are mangrove specialists, and three are introduced. Two families were represented by the anurans, five species from Eleutherodactylidae and one from Leptodactylidae. Parulidae was the most represented bird family, with five species, followed by Tyrannidae with three species. Columbidae, Cuculidae, Icteridae, Rallidae, and Vireonidae were all represented by two species. Ardeidae was also represented by at least two species, Egretta tricolor and likely a few herons and egrets (e.g., Butorides virescens, Bubulcus ibis, and Ardea herodias), whose similar vocalizations could not be differentiated by the authors. The remaining families were each represented by only one species and they are Accipitridae, Alcedinidae, Coerebidae, Mimidae, Pandionidae, Phasianidae, Picidae, Psittacidae, Scolopacidae, Strigidae, and Turdidae. Psittacidae was also likely represented by at least two species (e.g., Brotogeris versicolurus and Myiopsitta monachus) that could not be distinguished in the recordings. All the anuran species are invertivores and 14 (45%) of the birds as well. Of the remaining bird species, nine (29%) are omnivores, five (16%) are piscivores/carnivores, two (6.5%) are frugivores/nectarivores, and one (3%) is a granivore.

For thirty of these species, enough detections were identified to develop both pattern recognition templates and Random Forest models. The remaining six species (Actitis macularius, Buteo jamaicensis, Crotophaga ani, Gallus gallus, Mniotilta varia, and Setophaga adelaidae) were rarely detected or exhibited vocalizations that could not be easily modeled and are therefore only represented by validated pattern recognition. Eleutherodactylus coqui had the highest mean detection rate at 17% of nighttime recordings, followed by Coereba flaveola at 14% of daytime recordings. All other species were detected in less than 10% of either daytime or nighttime recordings. Parkesia noveboracensis was the most ubiquitous species and was detected at all 20 sites, followed by Coereba flaveola at 19, Melanerpes portoricensis at 18, and Ardeidae spp. and Quiscalus niger at 17. Eleutherodactylus brittoni and Setophaga adelaidae were the rarest species and were each detected at one site only.

Species accumulation curves show that the mean number of minutes required to reach 50, 75, and 100% of all the species detected at a site are 68, 308, and 1,396, respectively (Supplementary Figure 1A). But these times varied greatly among sites. Because all watersheds were not recorded for the same duration, the species accumulation curves were extrapolated to 10,000 mins for each watershed to predict a total of 35, 28, and 24 species in San Juan, Levittown and Ponce, respectively (Supplementary Figure 1B). Randomly sampling ten hectares between the three watersheds is expected to result in thirty to thirty-six detected species (Supplementary Figure 1C). Detections across all species and all sites dropped from a mean of 0.7 detections per recording to 0.3 detections per recording between the first deployment in March and the second in November. Further, there was a significant difference in the number of species detected at each site between the two sampling periods (t-test, mean difference = 3.9, p < 0.001). Of the thirty-six species detected between the two deployments, 34 were detected in March and only 29 in November. Of those that were not detected in November, Mniotilta varia, Setophaga americana, and Vireo altiloquus are migratory, while Setophaga adelaidae, Setophaga petechia, Margarops fuscatus, and Molothrus bonariensis are not. Buteo jamaicensis and Eleutherodactylus brittoni were detected only in November.

The remaining results apply only to the March 2017 survey, which occurred before hurricanes Irma and Maria. The top three linear models for all dependent richness variables are provided in Supplementary Table 3. Differences in avian richness across all sites were most strongly correlated with surrounding land cover and the minimum canopy height (Figure 2). Anuran richness increased logarithmically with the percent of surrounding land cover as mixed forest within 0.5 km. For birds, the total number of detected species at any site was most strongly and positively related to the minimum canopy height, which was also true for the number of invertivore birds. Endemic bird species increased linearly with the percent of mangrove cover within 5 km, while mangrove specialist species increased with surrounding freshwater forested wetland. Introduced species increased logarithmically with surrounding street density. Within other feeding guilds, carnivores/piscivores richness was negatively related to the percent estuarine scrub within 5 km, while frugivores/nectarivore richness was positively related to the percent developed open space within 0.5 km. Granivore richness was negatively related to surrounding urban cover and omnivore richness was positively and logarithmically correlated with surrounding mixed forest cover.

Figure 2. The strongest linear models to predict diversity among multiple groups of both bird and frog species in the urban mangroves of Puerto Rico. Land cover and canopy height were the strongest predictors, with urban variables most strongly predicting many of the bird groups.

None of the strongest linear models to predict richness metrics included the urban index and only one feeding guild, granivores, resulted in a significant difference in richness between the most urban and least urban sites at the island scale. Further, granivores are represented by only one species, Zenaida asiatica, which was absent in all five of the most urban sites and present at almost all the others. However, when comparing the local change in species richness between the minimum and maximum urban sites, six of the nine locations (BAH, LEV, MPN, PIN, PON, and TOR) exhibited higher bird richness at the maximum urbanness site (MAX), than at the minimum urbanness site (MIN). As a result, MAX sites have significantly more bird species than MIN sites (difference = 2.4, p = 0.01) (Figure 3), although the same comparison for frogs was insignificant (difference = 0.4, p = 0.3). Because all frog species are endemic and invertivores, the same result is true for frog endemics and frog invertivores. For birds, there were also more endemic (difference = 0.7, p = 0.02), introduced (difference = 0.8, p = 0.008), invertivore (difference = 1.2, p = 0.02), and omnivore (difference = 1.0, p = 0.01) species at the maximum urbanness sites in comparison to the minimum urban sites. Part of this is due to a significant dependence (Fisher’s test, p = 0.03) between the MAX urban sites and the presence of Psittacidae spp., which was present only at the MAX sites for four locations and was never present only at a MIN site. Although Psittacidae spp. was the only significant species with such a relationship, several other species showed similar, albeit, insignificant trends. Quiscalus niger and Tyrannus dominicensis were all present only at MAX sites for three locations, and never present only at a MIN site (Fisher’s test, p = 0.1). Five other species each were present only at a MAX site for two and one location(s), and never present only at a MIN site. On the other hand, Setophaga petechia and Rallus longirostris were only present at a MIN site for three locations and never present only at a MAX site (Fisher’s test, p = 0.1), while Crotophaga ani was only present at a MIN site for one location and never present only at a MAX site. Thus, for species that were present only at MIN or MAX sites and never at the other, thirteen were present only at a MAX site and three were present only at a MIN site, while seventeen were present at both MIN and MAX sites.

Figure 3. Differences in the number of species between the MAX urban and MIN urban sites at each of the nine locations. The most urban sites held significantly more species than the least urban sites.

Non-metric multidimensional scaling revealed some influence of urbanization on the bird and anuran community structure of the mangroves and this influence was again dependent on scale (Figure 4). At both the local and island-wide scales, the most urban faunal communities were more clustered than the least urban communities, suggesting the more urban faunal community structures are more similar to each other than are the less urban communities. Further, at the island-wide scale, ADONIS revealed a significant difference (p = 0.009) between the least urban and most urban communities, but the same test at the local scale was insignificant. The most urban sites were associated with the presence of C. flaveola, G. galeata, E. antillensis, E. coqui, and Q. niger, while M. antillarum, M. nudipes, P. haliaetus, and R. longirostris were associated with the least urban sites.

Figure 4. Non-metric multi-dimensional scaling of all sites based on the detection frequencies of their species, colored by island-wide relative urbanness (left) and local relative urbanness (right). Overlay vectors are land cover variables and species that were significantly correlated with and most responsible for the ordination. Sites are colored to the right as those with an urban index greater than the median (most urban) and those with an index value less than the median (least urban), while those to the left are colored based on their relative urbanness within each of the nine locations. Ellipses show 95% confidence areas of each group of sites. At both scales, the most urban sites are more clustered than the least urban sites, suggesting the former are more similar to each other, although analysis of variance in the ordination resulted in insignificant differences between urbanness groups at both scales.

Discussion

Previous studies on urban bird communities have consistently shown an influence of the urban landscape on faunal community structure and richness, which is largely attributed to the diversity, structure, and availability of surrounding habitats and resources. Studies in mangrove systems are much less abundant and show some consistency in those trends. This study, despite limitations related to unplanned for hurricane disturbances and acoustic identification, found both corroborating and contradictory findings in comparison to other studies. Primarily, that canopy height is most influential in explaining avian richness, but also that habitat cover and some urban variables play a role in structuring mangrove faunal communities, although scale is important in interpreting these results.

The mangrove avifauna richness of thirty-two species in the present study was five species more than that reported by Acevedo and Aide (2008), with fifteen of the same species shared between the two studies. Together, both studies suggest that the mangroves of Puerto Rico harbor at least forty-one bird species, thirty-two residents and nine migrants. Interestingly, Acevedo and Aide (2008) captured all nine migrants, while the present study captured only four. This highlights a disadvantage of acoustic avifaunal surveys, which are prone to miss non-vocal species, or species that mainly produce alarm calls, such as some non-breeding migrants. However, the present study detected fifteen additional resident species that were not detected by Acevedo and Aide (2008), many being relatively rare, cryptic, and/or potentially misanthropic and thus more likely detected by passive recorders than by visual point counts. Nonetheless, species-area curves suggest we identified nearly all of the vocal species inhabiting these forests during the observational periods.

Both the present study and Acevedo and Aide (2008) detected significant seasonal differences in the avifauna community, but while Acevedo and Aide (2008) detected higher richness attributed to migrants from August to February, the present study detected fewer species during the same period in November. Our lack of many of the migrants in November may be a result of timing, in which we missed the relatively short and fluctuating migratory period. However, only three of the eight species missing from the November survey are migratory, suggesting the absence of the other five is due to something other than normal seasonal changes in presence and/or vocalization behavior. Thus, it may be due to the forest disturbance associated with hurricanes Irma and Maria, which passed near and over Puerto Rico 2 months before the November survey, causing significant canopy loss (Branoff, 2020a). This is corroborated by findings from a more generalized (not mangrove specific) census of birds that resulted in significantly fewer birds counted after the hurricane (Wunderle, 2017). However, to further confound the patterns, four of the five non-migratory missing species from the present study are invertivores and the remainder is M. fuscatus, all of which were noted by Wunderle (2017) to increase after the hurricane. Thus, there is no clear explanation of the before and after hurricane differences and subsequent surveys are necessary to demonstrate the long-term resiliency of these communities to intense hurricane disturbances.

In simple linear regression, the richness of a few avian groups was most strongly explained by urbanization metrics. The percent of developed open space, defined as maintained green areas with some constructed materials but mostly low-lying vegetation, was the most important predictor for the richness of frugivores/nectarivores, while street density most strongly predicted the richness of introduced species and urban cover most strongly predicted the presence of the one granivore in the study. Both frugivores/nectarivores and introduced species increased with these urban metrics, which may be due to an abundance of food resources associated with urban environments. Other studies have noted the suitability of urban habitats for frugivores/nectarivores due to an abundance of both native and exotic flowers and fruits (Souza et al., 2019), while introduced species have also been found to be particularly abundant in streetscapes in Australia (White et al., 2005) and in Puerto Rico (Ruiz-Jaén and Aide, 2006; Suarez-Rubio and Thomlinson, 2009), again attributed to the presence of exotic vegetation from residential landscaping. In contrast, the one granivore in this study was only found in sites with less than 40% urban cover, despite other studies suggesting favorable conditions in urban habitats for these birds (Chace and Walsh, 2006).

At the local scale, pairwise comparisons of the most urban sites (MAX) with the least urban sites (MIN) resulted in significant differences in species numbers between the two. The maximum urban sites harbored around two more species per site than did the corresponding minimum urban sites for each location. This was true for six of the nine locations in the study. The same statistical differences were detected in the number of endemic, introduced, invertivore, and omnivore birds, each of which resulted in more species at the maximum urban site when compared to the minimum urban site. This may suggest that some species prefer the most urban site over the least urban site in each location, and the previously discussed patterns in omnivores and introduced species corroborate this. However, testing for dependence in presence/absence of individual species between MAX and MIN sites resulted in only one significant result. Psittacidae spp. was always present only at the MAX site when it was present at only one site at a location, which was four different locations. This is not surprising, as the two likely species that represent these detections, Brotogeris versicolurus and Myiopsitta monachus, have been noted to inhabit intense urban landscapes of Manaus, Brazil, and New York City (Bull, 1973; Souza, 2017). Twelve other species were also always present at the MAX site when they were present at only one site per location, compared to three species that were present only at the MIN site. Thus, when a gradient in urbanization is present within a mangrove forest, it seems more species are likely to be present in the most urban site than in the least urban site. The species which were most often present only at a MAX site include Quiscalus niger and Tyrannus dominicensis, which have previously been noted as successful urban inhabitants and should probably be considered synanthropes that benefit from these human landscapes (Ruiz-Jaén and Aide, 2006; Suarez-Rubio and Thomlinson, 2009). Q. niger for example is known to forage off of human garbage, while T. dominicensis utilizes urban street lights for easier insect capture (Raffaele et al., 2003). On the other hand, Setophaga petechia and Rallus longirostris were always only detected at a MIN site if it was detected at only one site per location, which agrees with previous studies that found these species to be urban-sensitive (Foin et al., 1997; Tremblay and St. Clair, 2011). Thus, these species should probably be considered misanthropes in urban mangroves.

Despite these clear patterns at both the island-wide and local scale, they often disagreed with one another between the scales. For example, there is no clear explanation for why endemics and invertivores were both more present in the MAX sites though neither were positively related to urban metrics at the island scale. One previous study has shown similar disjointed and inverse relationships between urban avifauna alpha and beta diversity at different scales (Norris, 2012). Thus, it is possible that multiple drivers are acting at different spatial scales, resulting in higher species richness in the most urban sites locally, but with little to no effect of urbanness on mangrove fauna richness across the island of Puerto Rico.

In this study, canopy height was often the most important covariate in explaining the variability of avian richness, and this covariate was also found to be important for Acevedo and Aide (2008), as well as for other mangroves (Goodenough, 2009) and for global forests in general (Feng et al., 2020). This pattern is attributed to greater micro-climate and niche space in taller canopies, although our analysis specifically tested for a relationship with the standard deviation of canopy height, which produced a relatively weak link. Other forest metrics, such as patch size and surrounding landscape matrix diversity have also been found to be important in urban mangrove bird communities (Mohd-Azlan and Lawes, 2011). Our results indicate no relationship between avian richness and surrounding vegetated habitat diversity but a positive relationship with the natural log of patch size (R² = 0.3, p < 0.05), which is the opposite of that found by Mohd-Azlan and Lawes (2011). This may suggest that these metrics are location specific or that different methodologies are utilized between the studies. In any case, maintaining high canopies in Puerto Rico’s urban mangroves may be the single most effective means to increasing avifauna richness, although this may be challenging in other areas where forests are still regenerating from previous development or where wood harvesting is still practiced (Mohamed et al., 2009; Branoff, 2018).

Mangrove specialists are of particular interest due to their dependence on mangrove habitats and thus their vulnerability to mangrove deforestation from urban development. Other studies have suggested mangrove specialists are capable of inhabiting urban mangroves (Ward, 1968; Noske, 1995, Lim and Sodhi, 2004), but these were nectarivores, which are absent in Puerto Rico’s mangrove specialists. The only two mangrove specialists present in the most urban mangroves of this study were Zenaida aurita, an omnivore, and Parkesia noveboracensis, an invertivore. Both species are defined as “mangrove dependent” rather than “mangrove restricted”, suggesting they can utilize other habitats in addition to mangroves, which may explain their presence in Puerto Rico’s urban forests. Further, Z. aurita is an omnivore and may benefit from the previously mentioned variety of foraging opportunities offered in urban mangroves, making these forests just as suitable, if not more so, than less urban mangroves. On the other hand, P. noveboracensis is an invertivore and was ubiquitous across all the surveyed sites, suggesting it is relatively flexible in habitat preferences, regardless of urban intensity. The two mangrove restricted specialists, Rallus longirostris and Setophaga petechia, were absent from the most urban sites, suggesting they are relatively more misanthropic and adverse to inhabiting the urban mangroves of Puerto Rico. Additionally, endemic bird richness increased with surrounding mangrove cover, suggesting mangrove habitat is particularly important for these species, despite none of the endemic birds being mangrove specialists. This further highlights the importance of mangrove cover as habitat for a wide range of avifauna and particular attention and additional conservation resources may be warranted for these species where urbanization threatens mangrove deforestation.

Anuran richness could not be analyzed among multiple groups or dietary guilds, as all five species are endemic invertivores. Thus, the only predictor for any of the anuran groups was the percent cover of surrounding mixed forests, which was positively correlated with richness. As with birds, there was little evidence for an island-wide influence of urbanization on anuran richness. This is not surprising, as previous studies have also shown no effect of traffic noise on anuran acoustic activity in Puerto Rico (Herrera-Montes and Aide, 2011), despite declining numbers in lowland forests on the island (Campos-Cerqueira and Aide, 2017). Interestingly, none of the anuran species in this study have previously been described in mangrove habitat and are typically associated with palustrine, not estuarine, wetlands. This is likely because most anurans are poor osmoregulators, especially as larva (Gomez-Mestre et al., 2004). However, the direct development of Eleutherodactylus spp. probably negates poor osmoregulation, allowing them to maintain their lifecycle despite the elevated salinity in mangrove systems and with no effect from surrounding urbanization. As the Caribbean has been identified as one of five biodiversity hotspots predicted to undergo relatively rapid urbanization in the twenty-first century (Seto et al., 2012), mangrove avian and anuran faunal communities there may experience shifts toward less richness and more homogeneity. This study showed some distinction of faunal communities between the most and least urban mangroves of Puerto Rico, although scale played an important role. Island-wide, faunal community structure was distinct between the least urban and most urban sites but in avian richness, it appears urbanization is not as important as forest structure and minimum canopy height. At the local scale, however, there were no differences in faunal community structure between maximum and minimum urbanness, but the maximum urbanness sites held significantly more bird species than their counterparts. Although the extra species at the maximum sites were not always consistent across locations, there is some indication that numerous avian omnivores and synanthropes explain the elevated richness of the most urban mangroves. On the other hand, some misanthropes, as well as endemics and mangrove specialists in general will likely benefit from minimal urban influence and more expansive mangrove forests. At the same time, highly urban forests can be managed to optimize cultural ecosystem services that take advantage of both avian and anuran use of these habitats. This includes recreational services such as bird watching, and educational services focused on ecology. Both strategies will be important for the long-term sustainability of these social-ecological systems.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://osf.io/uk2fd/.

Author Contributions

BB contributed to the study conception, design, implementation, completion of statistical analysis, also performed analysis and acoustic pattern species identification, and wrote the manuscript. MC-C performed training, verification, validation, and revision of the manuscript. Both authors contributed to the article and approved the submitted version.

Funding

Portions of this work were funded by the United States Department of Agriculture’s Forest Service International Institute of Tropical Forestry.

Disclaimer

The views expressed in this manuscript are those of the authors and do not necessarily represent the views or policies of the United States Environmental Protection Agency. Any mention of trade names, products, or services does not imply an endorsement by the United States Government or the United States Environmental Protection Agency. The EPA does not endorse any commercial products, services, or enterprises.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Orlando Acevedo-Charry for his assistance in the field and T. Mitchell Aide for his assistance in acoustic analyses.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.670288/full#supplementary-material

Footnotes

1. ^ https://arbimon.rfcx.org/project/mangroves-of-puerto-rico/dashboard

References