The Structure and Composition of Puerto Rico’s Urban Mangroves

This study characterizes the structure and composition of mangrove forests across urban gradients in Puerto Rico. It then uses a suite of hydrologic, water chemistry, and land cover variables to test for the relative importance of urban intensity alongside flooding and water chemistry in explaining observed variability in forest structure and composition. Three separate statistical tests suggest a significant but limited influence of urbanness on forest composition and structure. In the most urban sites, the diameters of the largest trees were 27% larger, but all structural measurements were best explained by surface water chemistry, primarily nitrogen concentrations. Concentrations of ammonium and total Kjeldahl nitrogen best explained stem density, tree girth and canopy height. The most urban forests also contained 5.0 more species per hectare, on average, than the least urban forests, and simple regression suggests that urban metrics were the most powerful predictors of forest composition. The most urban forests were more dominated by Laguncularia racemosa, while both Avicennia germinans and Rhizophora mangle were found to be less abundant in the most urban sites, a trend that may be linked to the influence of precipitation and tidal connectivity on porewater salinity across the urban gradient. In multiple regression, no statistical difference was detected in the importance of surrounding land cover, flooding, or water quality in explaining the variance in either composition or structural metrics. This suggests that while a given forest metric may be strongly linked to either land cover, water quality, or flooding, all three are likely important and should be considered when characterizing these forests. With more human dependents in urban areas, the provisioning of important ecosystem services may be influenced by land use variables in addition to the more commonly measured metrics of water chemistry and flooding.


Introduction
Urbanization has been an important contributor to forest disturbance at the turn of the twentieth century, with trends indicating an increasingly important role in tropical coastal areas over the next few decades [1][2][3][4]. The conversion of forests to developed lands usually follows economic transitions that favor industrialized societies and continuing transitions may in some cases support regrowth, albeit with novel forests [5,6]. As a result, forests in urban and surrounding lands are often represented by unique anthropogenic influences and ecological traits [7][8][9].
Characterizations of urban forests often consist of lower stem densities and larger individuals, with stands exhibiting increased edge openness and regeneration failure [8]. These forests are further generalized as having higher floral diversity, usually due to nonnative species from residential gardens and municipal landscaping. As a result, the function of these systems is also novel, as evidenced by altered community composition, biogeochemistry, productivity, and resiliency [10,11]. This, in turn, gives way to an adjusted provisioning and valuation of ecosystem services [12,13]. Much of this is known from the study of terrestrial systems, with comparatively little understood of mangroves, the dominant forested biome of tropical, sub-tropical, and warm temperate estuaries.
Globally, mangrove coverage in the largest cities is decreasing faster than overall rates from corresponding countries [14]. In some cases, this loss results in fragmented forests consisting of novel species assemblages and size classes [15][16][17][18]. There are cases of expanding urban mangrove coverage [19,20], although the young forests may contain less biomass [21]. These studies suggest a systemic influence of urban land-use on mangrove structure and composition, but most were done independently of hydrology or water chemistry, which have long been recognized as powerful influences on the same structural and compositional metrics [21][22][23].
Flooding metrics of hydroperiod and flood frequency are important for mangrove seedling survival and adult growth [24,25]. Water salinity, pH, dissolved oxygen, phosphorus, and nitrogen concentrations are also known to influence mangrove physiology and forest structure and composition [26][27][28][29][30]. While there are some anecdotal connections between mangrove hydrology, water chemistry, and urbanness [31][32][33], quantified connections between these components are absent or weak, leaving little definitive evidence for the hypothesis that urban land use influences mangrove forests through hydrology and water chemistry.
In Puerto Rico, Brandeis,et al. [34] provide a limited structural inventory of the mangrove and non-mangrove forests of the island's largest city, San Juan, based on a small sample area, but there are no conclusions regarding the influence of urban areas on this structure, and much less is known of the urban mangroves outside of the San Juan metropolitan area. Overall, urban mangrove coverage in Puerto Rico has largely failed to expand during the last quarter of the twentieth century, despite a 12% increase in mangroves across the island [35]. Furthermore, urban forests were found to have fewer and smaller forest fragments compared to more rural sites. However, while there is some limited information on the structure and composition of these forests, there is little empirical evidence tying any observed changes in mangrove ecology to specific components of urban landscapes. Furthermore, as the 2017 hurricane season significantly disturbed these forests [36], it is important to have a predisturbance benchmark by which to monitor and compare subsequent recovery among the urban and non-urban forests.
In addition to the above shortcoming, there are no universally accepted definitions for "city" or "urban", and thus little consensus on what defines an urban forest and our understanding of their structure and function [37,38]. Many of the above described studies use qualitative comparisons of urban vs. non-urban but provide no definitions for these characterizations, and a very limited few attempts to quantitatively define urbanness through landcover (e.g., impervious surfaces) or demographics (e.g., population density). As a result, it is difficult to repeat or synthesize studies to understand systemic patterns in urban forest ecology. Thus, urban ecologists increasingly recognize the need for quantitative metrics of urbanness that can be represented as an urban gradient [11,39,40]. Such metrics not only allow for a more causal understanding of urban ecosystem patterns towards hypothesis testing, but also for more repeatable studies across multiple urban landscapes.
This study uses a combination of ground-based and remote sensing measurements to address the following objectives: 1) Characterize the structure and composition of mangrove forests along quantified urban gradients in three watersheds of Puerto Rico, 2) Test the hypothesis that urban forests are characterized by lower stem densities and stand biomass, larger individuals, and higher floral diversity when compared to less-urban counterparts, and 3) Test for the combined importance of urbanness, hydrology, and water chemistry on mangrove forest structure and composition through multiple regression. Addressing these objectives will be an important component of mangrove management under increasing pressure from urbanization in the Caribbean and worldwide [4,[41][42][43].

Study Location
We focused on twenty mangrove sites (one-hectare each) from Branoff [45] and distributed in the coastal areas of Puerto Rico. These sites were established according to a previously developed urban index [44], such that they fell within the greatest range of urbanness ( Figure 1, Table 1). The urban index is a unitless metric of urbanness that incorporates information on urban land use, vegetation and open water coverage, population density, and road length, where a value of 0 represents the least urban site and a value of 100 is the most urban site (Table 1). Fourteen sites fall within the San Juan Bay Estuary and three sites each Branoff and Martinuzzi Page 3 Forests. Author manuscript; available in PMC 2021 October 21. fall within watersheds in Ponce and Levittown ( Figure 1). Sites in the San Juan Bay Estuary are named for their waterbodies (BAH is the San Juan Bay, MPN and MPD are the undredged 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). In Ponce and Levittown, mangrove area is comparatively smaller and so sites are named for their watershed only: LEV is Levittown and PON is Ponce. To ensure we sampled the greatest feasible range of urbanness in a given area, two sites per waterbody in the San Juan Bay Estuary and three sites each in Ponce and Levittown were selected, each representing the relative minimum, mid, or maximum urban index values for that area. In the San Juan Bay Estuary, two sites were located in the minimum and maximum levels of urbanness in each waterbody, as denoted by "MIN" and "MAX" following waterbody abbreviations in site names. In Levittown and Ponce, three sites were placed at the minimum, median, and maximum urbanness levels of each watershed, as denoted by "MIN", "MID" and "MAX" following watershed ("PON" and "LEV") abbreviations. Thus, "MIN", "MID" and "MAX" in site names are relative to local urbanness levels only, not island wide. All sites are bound on one side by a waterbody and most sites are roughly 100 m by 100 m, with some asymmetry due to non-linear coastlines. In some cases (BAHMIN, BAHMAX, MPDMIN, MPNMIN, MPNMAX, TORMAX) forests were constricted by natural or anthropogenic features and these sites were extended along the shoreline to compensate for a lack of forest less than 100 m from the shore.
Human use of the mangroves is minimal at most sites, although harvesting of the crustacean, Cardisoma guanhumi, is sparsely evident at some and forest use for grazing and corralling of horses, pigs, and dogs is also apparent at MPNMAX and LEVMAX. There was no evidence of wood harvesting at any of the sites. For more information on the urban index and plot information see Branoff [44].
The mangroves of San Juan were most recently described by Brandeis,Escobedo,Staudhammer,Nowak and Zipperer [34] as being dominated by Rhizophora mangle, although the authors point out this is likely due to sampling methods and that Laguncularia racemosa best characterize the forests and represents more biomass than any other species.
Very little has been published on the mangroves of Levittown, which surround an estuary composed primarily of an artificial tidal lagoon constructed to drain surrounding settlements and connected to the ocean through a tidal creek and permanent inlet. Unlike the other two sites, the mangroves at Ponce are largely unconnected and do not share the same estuarine conditions among them. Additionally, Ponce receives a much drier southern climate in comparison to the northern sites, with a median annual rainfall of 755 mm versus 1600 mm on the northeastern coast [44]. Like Levittown, very little is known about the mangroves in Ponce, although one study did characterize the hydrology of Punta Cabullones [45].
All sites were classified as one of four hydro-geomorphologies following aerial imagery, previous studies, and the results of water level monitoring from Branoff [44]. These classifications are embayments (i.e., ocean, bay, or lagoon) or canals, and either open or partially restricted to tidal exchange. Due to a lack of tidal connectivity, water levels in partially restricted sites respond more strongly to rainfall events in comparison to open sites [44] and this is the primary criteria for their hydro-geomorphic classifications although the Branoff  following diurnal water level ranges also reflect respective tidal connectivities. Ocean and bay sites always have direct tidal exchange and include the two San Juan Bay sites BAHMIN and BAHMAX, and the ocean site in Ponce, PONMIN. Hydrographs of these sites show mean diurnal ranges of 40, 47, and 14 cm, respectively [44], reflecting that of offshore buoys in the Atlantic Ocean (48 cm) and Caribbean Sea (20 cm) [46,47]. Lagoon sites, however, may include embayments that are partially restricted to tidal flow. The only open lagoon sites are in Torrecilla lagoon (TORMAX and TORMIN), due to its dredged mouth at Boca de Cangrejos (Ellis 1976), as evident in a mean diurnal range of 41 cm [44]. Closed lagoons in the San Juan Bay Estuary include Piñones (PINMAX and PINMIN) and San José (SANMIN and SANMAX), whose mean diurnal ranges are reduced to 3 and 4 cm, respectively [44]. Elsewhere, closed lagoons include the Levittown lakes (LEVMAX) in Levittown, and in Ponce the Salinas lagoon (PONMID) and the Port of the Americas (PONMAX), with mean diurnal ranges of 0, 5, and 3 cm, respectively [44]. Canal sites may also have open or partial tidal exchange. Open canals are the east end of Suárez canal, SUAMAX, which shares the dredged connection with Torrecilla and has a mean diurnal range of 17 cm, the west and dredged portion of the Caño Martin Peña, MPDMIN and MPDMAX, both with mean diurnal ranges of 47 cm, and the mouth of the Río Cocal in Levittown, LEVMIN, with a mean diurnal range of 13 cm [44]. Partially restricted canals are that connecting the Río Cocal to the Levittown lakes, LEVMID, with a mean diurnal range of 12 cm, the undredged eastern portion of the Caño Martín Peña, MPNMAX and MPNMIN, each with diurnal ranges of 27 cm and the portion of the Suárez canal restricted by the Baldorioty expressway (SUAMIN), which shares the San José range of 4 cm [44].

Ground-Based Measurements
Ten subplots within each one-hectare site were sampled, each being a 5 m radius circle of 78.5 m 2 . A 5 m radius plot design was used as recommended for small stem forests [48]. Plots were established by fixing the end of a 5 m rope at ten randomly generated coordinates within each site and extending the rope until tight. These center points were located using a Garmin eTrex 10 global positioning system with an accuracy of ± 3 m. All woody plants greater than 1 cm in diameter at the breast and within the 5 m radius were identified and their diameter measured with a diameter measuring tape at a height of 1.4 m, or diameter at breast height (dbh). For individuals of R. mangle whose primary prop roots joined the trunk at heights greater than 1.4 m, trunk diameter was measured just above the confluence of the prop roots where a true mainstem existed. Branches from mainstems of all species were also measured if they originated at heights less than 1.4 m and if their dbh also exceeded 1 cm. Finally, in addition to the surface water chemistry metrics described below, we also took porewater salinity measurements at each site by extracting a 5 mL sample of porewatersediment mix through the use of a small diameter tube connected to a 50 mL syringe similar to that described by McKee, et al. [49]. This was performed at depths of 0, 10 and 20 cm, and the resulting extraction was allowed to settle until a clear sample of water could be measured for salinity, in ppt, via a VEE GEE STX-3 refractometer. In some cases, a clear sample of water could not be obtained from the sediment and no salinity was recorder for that depth at that site. Branoff  Aboveground biomass estimations were calculated for each tree through allometric equations using the measured dbh and wood specific gravities when available. Halophytic plant types and salinity tolerances, when available, were provided by Santos, et al. [50]. For the three true mangrove species (A. germinans, L. racemosa and R. mangle), equations were species and dbh specific as derived from three separate sources on Caribbean mangroves [51][52][53]. The mean of these three values was used and when no value was available for greater size classes, a general equation for mangrove habitats was used from Chave, et al. [54]. This equation was also used for non-mangrove species in combination with specific gravities derived from Reyes, et al. [55]. All ground measurements used to characterize mangrove structure and function are described in Table 2.

LiDAR
Airborne LiDAR data were collected over parts of Puerto Rico in March of 2017 using the NASA G-LiHT (Goddard's LiDAR, Hyperspectral and Thermal) imaging system [56], as part of the US Department of Energy, Next-Generation Ecosystem Experiments-Tropics project (https://ngee-tropics.lbl.gov) ( Figure 1). LiDAR data provide detailed information on forest 3D structure and therefore can be used to quantify forest structure across urban gradients [57,58]. Among our sites, LiDAR data from the GLiHT campaign were only available for the San Juan Bay Estuary sites and the LEVMID site. Data were collected at a flying altitude of 335 m above ground and ~12 pulses m −2 . Point clouds were processed in R version 3.6.1 [59] and the lidR [60] and rLiDAR [61] packages. The LiDAR point clouds were clipped to the 5 m radius plot boundaries using the lasclip function and height and canopy metrics extracted from the LASmetrics function. Common LiDAR metrics returned by LASmetrics were derived for each plot, including: mean height, standard deviation of height, height percentiles, skeweness, canopy cover, and canopy density (see Table 2 for more information). We used the 90th height percentile to represent the maximum canopy height throughout the analysis to avoid outliers associated with actual maximum return heights [62]. The geolocation error of the field plots (± 3 m) should not be a problem, as the LiDAR data was used to characterize the sites, and not to relate the LiDAR data and field data at the plot level.

Land Use, Flooding, and Water Chemistry
To explain the variability in the above described forest structure and composition metrics, we used several land use, flooding, and water chemistry variables as outlined in Table 3. Because analysis at multiple scales is important and ecological effects of land cover are typically most pronounced within 1 km [63,64], land cover data were cropped, masked, and sampled using circles of radius 50-1000 m around each plot. For each circle we calculated the extent of three classes: urban (including the classes Impervious and Developed Open Space, respectively), green and blue area (represented by the sum of all vegetation and open water classes), and mangrove, corresponding to estuarine scrub/shrub wetland and estuarine forested wetlands. This was accomplished through the st_buffer, st_crop, and st_mask functions of the sf package [65] and the getValues function from raster [66].
Road length was calculated by clipping a roads layer and summing the length of all resulting road segments. Population density was calculated from the 2010 Census and the land cover Branoff Forests. Author manuscript; available in PMC 2021 October 21. and road layers previously described by assuming people lived only in non-road impervious surfaces. These variables describing land cover, population, and roads are included in a single urban index that represents the relative urbanness of each site [44]. The urban index first standardizes each variable from 0 to 100, 0 being the least urban and 100 being the most urban value, then computes an average of each normalized variable to give the single urban index representing a combination of all variables [45]. In addition, we used information on site hydrology and water quality derived from water level models and surface water quality measurements representing the period from 2012 to 2017 at the same study sites [44]. Water level models were constructed from water level observations recorded by data loggers at each site and precipitation observations at nearby weather stations. Tidal harmonics models were used to extract tidal constituents and moving sums were used to model precipitation contributions to observed water levels. These were then used to reconstruct water levels over a five-year period in which flooding conditions were assessed every hour by comparing the predicted water level to the elevation of the forest as determined from digital elevation models [45]. Hydrology metrics include average depth (m), proportion of time flooded (fraction), mean daily flood frequency (floods day −1 ), mean flood length (days), number of days with at least one flood per year (days), flooded hours per year (hours) and dry hours per year (dry). Water chemistry measurements were made in surface waters within 1 km of each site on a monthly and bi-annual basis over the same five-year period by the San Juan Bay Estuary Program and are available at https://estuario.org/. These measurements were only available for the fourteen sites in San Juan. Monthly measurements are obtained in-situ by a handheld sonde and are dissolved oxygen (mg/L), pH, salinity (PSS), specific conductivity (ms/cm), and temperature (°C). Bi-annual measurements are obtained by water samples from surface waters and laboratory chemical analyses. These include ammonium (mg/L), total Kjeldahl nitrogen (mg/L), nitrate and nitrite (mg/L), and phosphorus (mg/L). A detailed description of the water quality and hydrology measurements is given in Branoff [44] and scripts and results for this analysis can be found at https://github.com/BBranoff/Urban-Mangrove-Hydrology.

Analyses and Statistical Testing
All numerical analyses and statistical tests were performed in R version 3.6.1 [59]. Except for regression analyses, statistical tests were performed at the site and watershed levels using values from each plot (n = 200; 10 plots for each of the 20 sites). Comparisons across sites and watersheds were done through an ANOVA as calculated by the aov function of the stats package [59], and pairwise tests were performed through a Tukey Honest Significant Differences test as computed through the TukeyHSD function in the same package. Nonmetric multidimensional scaling (NMDS) was performed to observe the relative differences in forest structure and composition between sites and hydro-geomorphologies (i.e., open embayment, restricted embayment, open canal, restricted canal), and to observe the potential contribution of the urban index to any separation and/or grouping of forests. This was performed twice, once on compositional measurements (e.g., species diversity, percent of biomass as R. mangle, etc.) and again on structural measurements (e.g., maximum dbh, stem density, biomass etc.). Functions for this analysis are from the vegan package [70]. NMDS was first performed through the metaMDS function with the default of two dimensions. Results were then rotated through the MDSrotate function so that the horizontal axis of the NMDS was aligned with the vector describing the urban index. Vectors describing the various measurements were then calculated through the envfit function, which finds the projection of the maximum correlation between each variable and the ordination. Finally, ellipses describing the 95% standard error confidence area for each hydro-geomorphology were drawn using the ordiellipse function. Sites in the NMDS were grouped into hydrogeomorphic characterizations as reported in Branoff [44] to account for any potential groupings due to differences in flooding dynamics.
The contribution of environmental predictor variables (e.g., urban cover, proportion of time flooded, ammonium concentration, etc.) on forest structure and composition (e.g., stem density, biomass, number of species, etc.) was analyzed through both simple and multiple regression. To single out the most important environmental influences on forest metrics, simple regression models were constructed through the lm function [59] in the form y~x, y~ln(x), and y~ln(x+1) when environmental predictor variables included values less than or equal to 0. The highest performing models were selected as those with the highest R 2 value whose p-value was lower than 0.05. Regression was done using mean site values to better highlight general trends across sites. These are accompanied by boxplots comparing the same forest composition and structure metrics between the least urban, urban, and most urban plots.
To examine the combined and relative importance of urban, hydrology, and water chemistry variables in explaining variability in forest structure and composition, multiple regression models were constructed by including one variable from each of the three potential influences of land cover, hydrology, and water chemistry. These models took the form of where Y is a given response variable of forest structure and composition, a is the intercept of the model, b are the predictor specific slopes, X are the predictors, and the subscripts 1, 2, and 3 represent groups of predictors from land cover, hydrology, and water chemistry. Models were constructed in the same way as simple regression, with both raw predictor variables as well as their natural logarithm to test for both linear and non-linear relationships. Bayesian Information Criteria (BIC) were calculated for all candidate models of the same response variable through the glmulti function of the same package name [71].
BIC was used because it imposes a higher penalty than AIC for models with multiple predictors, thus selecting for the simplest model. Collinearity in regression predictors was tested for through the ols_vif_tol function of the olsrr package [72], with any variance inflation factor (VIF) values greater than 4 prompting closer inspection of the model. The top performing significant model (p < 0.05) for each response variable were selected based on the lowest BIC value. The relative importance of each variable to final models, as determined by its contribution to the overall R 2 , was calculated through the calc.relimp function of the relaimpo package [73].
Graphs were produced through ggplot and the geom_violin, geom_bar, and geom_point functions for violin, bar and scatter plots, respectively [74]. Violin plots represent data at the plot level. Scatter plots represent data at the site level, and simple regression models were Branoff  plotted through the stat_smooth function using the formula y~x, unless it was outperformed by the formula y~ln(x) or y~ln(x+1) as determined by the highest R 2 value. R scripts as well as all required raw data can be found on the Open Science Framework project page: https:// osf.io/uk2fd/. This includes individual tree measurements and plot level characterization of land cover, water quality, and flooding metrics, as well as LiDAR metrics at the plot level. The 2017 LiDAR data for Puerto Rico are available at https://gliht.gsfc.nasa.gov/.

Structure and Composition of Mangroves
Nine-thousand three-hundred and eighty-one stems belonging to 7250 mainstems were identified to a total of 30 different species, resulting in an overall stem density of 5997 ± 291 stems per hectare ( urban site, was distinct in harboring high mangrove diversity at each plot, with an average of 2.7 species and was significantly more diverse than seven (35%) of the other sites (ANOVA; p < 0.05). LiDAR data were only available for sites in the San Juan Bay Estuary and for LEVMID (Figure 3

Non-Metric Multidimensional Scaling
The NMDS showed little separation of sites or hydro geomorphologies based on the compositional and structural metrics, suggesting all sites are relatively similar in composition and structure ( Figure 4). Stress values in NMDS scaling were 0.104 for compositional measurements and 0.13 for structural measurements, suggesting fair representations of actual data [75]. Although all variables were statistically significant when regressed on the ordination axes, the coefficient of determination (R 2 ) of these fits were lowest for the urban index (R 2 = 0.03 for composition and 0.05 for structure). Species diversity, percent stand biomass composition of L. racemosa, height, dbh, and canopy cover metrics were most aligned with the urban index, suggesting a correlation between these variables.

Simple Regression
In simple regression, forest composition was best described by models including the urban index, population density, and porewater salinity, whereas structural metrics were best explained by models involving flooding and surface water chemistry predictors ( Figure 5). The percent stand biomass composition of A. germinans decreased with increasing urban index, while that of L. racemosa and R. mangle decreased and increased, respectively, with increasing porewater salinity. Furthermore, while the percent stand biomass composition of L. racemosa was significantly greater in the most urban sites than the least urban sites, that of R. mangle showed the opposite trend and was greatest in the least urban sites. Moreover, the overall tree species diversity increased, while mangrove diversity decreased with increasing surrounding population density within 200 m. The significance of this model, however, was due to the relatively high species numbers at MPNMAX and could not be repeated when this site was removed from the test. Branoff Forests. Author manuscript; available in PMC 2021 October 21.
The other response variables, mostly structural, were best explained by water nitrogen and phosphorus concentrations ( Figure 5). Mean height, the 90th height percentile, the standard deviation of the height, and the maximum dbh all increased with increasing concentration of total Kjeldahl nitrogen in surface waters. The height skewness decreased with increasing surface water phosphorus, and the stem density decreased with increasing surface water ammonium. Both basal area and biomass decreased with higher surface water temperature. Finally, canopy cover increased with increasing surrounding water coverage within 200 m, and canopy density increased with higher daily flood frequency.
The standard deviation of height was significantly higher in the most urban sites when compared to the least urban sites (t-test, difference = 0.6 m, p < 0.01), as was the maximum tree diameter (t-test, difference = 5.7 cm, p < 0.01). Additionally, the most urban sites held more species than both the urban (t-test, difference = 0.7 species per plot, p < 0.01) and least urban sites (t-test, difference = 1.3 species per plot, p < 0.001).

Multiple Regression
To test for the influence of multiple variables on forest structure and composition, fifty-two thousand candidate models were constructed from all combinations of one, two, or three variables, one each from the groups of land cover, hydrology, and water chemistry metrics (Table 5). Top ranked models by BIC almost always included all three variables and resulting R 2 values averaged 0.74. Although the mean importance of land cover variables in explaining variance in forest composition metrics was highest at 54%, it was not significantly higher than water quality (mean = 32%) or flooding (mean = 26%) (ANOVA, p > 0.25). There was also no significant difference in the relative importance of land cover (mean = 33%), flooding (mean = 19%), or water quality (mean = 32%) variables in explaining structural metrics (ANOVA, p > 0.2). Specific urban variables of population density, road and highway density, and urban cover held an average relative importance of 28%, compared to 40% for water chemistry and 32% for flooding, again with no statistical difference between the three.

Discussion
Results from three separate statistical methods (i.e., Non-metric multidimensional scaling, simple regression, and multiple regression) suggest an indirect influence of urban land use on mangrove forest composition and structure in Puerto Rico, and that this influence is shared between flooding dynamics and water chemistry. Non-metric multidimensional scaling showed a weak influence of the urban index in the ordination of sites based on both compositional and structural metrics. It also showed a weak separation of individual sites and limited grouping, suggesting sites are relatively homogenous in structure and composition and not distinguished by their locations or hydro-geomorphologies. Still, simple regression models showed that forest composition was mostly explained by urban variables (i.e., population density and the urban index) and demonstrated increasing tree diversity but decreasing mangrove diversity with greater urbanness. In structural metrics, however, surface water nitrogen and phosphorus concentrations were the most common top ranked predictors and their shared influence on mangrove forests was further evident in Branoff Forests. Author manuscript; available in PMC 2021 October 21. multiple regression. This suggests that urban lands may be influencing these forests, albeit indirectly through pathways that are relatively weakly linked to the urban index.
Tree size, as indicated by the maximum and mean diameters as well as the 90th height percentile and mean height at each site, were most strongly predicted by surface water nitrogen in both simple ( Figure 5) and multiple regression (Table 5). In the case of diameter, this translates to trees with higher maximum stem diameters in the most urban forests, which were 27% larger than in the least urban forests. This may be because the most urban mangroves are permanently or periodically flooded (0.001-2.5 times per day) by urban waters, which likely carry elevated nitrogen and phosphorus concentrations that have been tied to roads and sewage [76][77][78]. Such nutrient enrichment has been shown to increase girth but decrease height of mangrove seedlings [79], while the height of dwarf mangroves was found to increase with both nitrogen and phosphorus fertilization [80]. Thus, these nutrients likely influence urban mangrove growth depending upon antecedent conditions and their interaction with other water quality and flooding conditions. Although a previous study at these sites found no linear link between the urban index and surface water nitrogen and phosphorus, it did suggest that urban wastewater effluents may be contributing to elevated concentrations of these elements in the most urban mangrove waters, which held three times as much nitrate and nitrite as the least urban waters [44]. Thus, while there is a relatively weak direct relationship between the urban index and mangrove size, scattered point source wastewater effluents in the most urban sites may indirectly contribute to overall larger trees there.
Despite the positive influence of nitrogen on individual tree size, overall forest structure as measured by stem density was found to be negatively related to ammonium in both simple ( Figure 5) and multiple regression (Table 5). This concomitant reduction in stem density with increasing tree size has also been observed in other forests of the Caribbean [51,[81][82][83], and can be expected from a self-thinning dynamic in which larger trees outcompete smaller surrounding trees, thus reducing stem density and increasing average tree size as the forest matures [51]. Yet, this theory would also predict an overall increase in forest basal area and biomass as stem density is reduced, a pattern we could not detect in the present study (Biomass ~ Stem Density: p > 0.1, R 2 < 0.05; Biomass ~ Ammonium: p > 0.3, R 2 < 0.1). One study in South Florida suggests that this tradeoff between stem density and biomass occurs only in the absence of stress [82]. Thus, the absence of such a pattern in the present study may be indicative of stress that is limiting forest basal area and biomass, despite sufficient nutrient conditions and larger trees. In fact, water temperature was the strongest predictor of both basal area and biomass in simple and multiple regression and physiological studies have shown an inhibitive effect of high temperatures on mangrove productivity, especially above 25 °C [84][85][86]. Thus, although the largest trees are larger in urban nitrogen enriched forests, it may be that while overall stem density is driven down by this dynamic, basal area and biomass are driven down by elevated water temperatures (> 30 °C). Perhaps as a result of these confounding factors, there are no differences in forest biomass between the most urban and least urban forests of Puerto Rico.
Compositional metrics were more consistently tied to urbanness than were structural metrics. Overall tree diversity was positively correlated with population density (Figure 5), Branoff  which agrees with terrestrial urban forests [8]. In the present study, this seems to be driven by the presence of non-mangrove, and in some cases non-native species, found only in the most urban sites, including Andira inermis, Mammea americana, Paullinia pinnata, Pavonia fruticosa and Tabebuia heterophylla. The presence of these species within the most urban mangrove forests may be a result of a combination of their increased use in surrounding municipal and residential landscaping [87] as well as lower salinity and inundation due to restricted tidal influence in highly engineered urban landscapes [44]. Little is known on the effect of urban land use on mangrove diversity, but other studies have pointed out the tendency of L. racemosa to form monoculture stands in post-disturbance shorelines [15,88]. This agrees with our results, in which the percent stand biomass composition of L. racemosa increased with population density ( Figure 5), resulting in roughly 40% more biomass represented by this species at the most urban sites compared to the least urban sites. Contrastingly, the contribution to mangrove biomass by A. germinans was negatively related to the urban index and that by R. mangle was significantly lower in the least urban sites compared to the most urban sites. For both L. racemosa and R. mangle, this may be tied to salinity, as both also showed a strong trend with porewater salinity, suggesting that changes in urban hydrology may be driving the observed compositional metrics.
Previous studies have largely hypothesized that observed patterns in urban mangrove species composition are primarily due to changes in hydrology from the urban environment [31,33]. Indeed tidal influence and connectivity are critical to mangrove zonation and competition with non-halophytes [89,90]. Thus, it was hypothesized in the present study that changes in the mangrove community along the urban gradient could be tied to changes in tidal connectivity and hydrology. A previous study at these sites did find evidence of reduced tidal amplitudes along the urban gradient as well as immediately upstream of canal restrictions or dredging projects [44]. The present study, however, found mixed evidence in support of this hypothesis. While mangrove composition did change along the urban gradient and a connection to porewater salinity was also observed, these changes were not consistently tied to changes in hydrology or to infrastructure projects. For instance, there was no significant difference in the number of species per plot or the percent stand biomass composition of A. germinans or L. racemosa between the dredged portion of the Caño Martín Peña at MPDMAX and the un-dredged portion roughly 500 m away at MPNMIN, which exhibits a near 50% (20 cm) reduction in the diurnal tide range. The same is true for SUAMAX and SUAMIN, the two sites on either side of the Suárez canal restriction at the Baldorioty expressway, which also exhibit a 14 cm difference in diurnal tide ranges between the two. Thus, changes in hydrology as a result of infrastructural projects alone cannot account for the observed change in compositional metrics along the urban gradient. Furthermore, while inundation parameters and soil porewater salinity were important in describing variability in forest composition under multiple regression, this importance was shared with and often shadowed by that attributed to surrounding urban metrics. As a result, it cannot be definitively concluded that hydrology or infrastructure projects are the most important drivers of urban mangrove floral community composition. Instead, there may be an alternative selection pressure in urban environments (e.g., residential and municipal landscaping) that drives the observed patterns, which may have implications for the Branoff  suitability of urban mangroves to provide habitat to a wide variety of taxa that often seek specific floral species [91].
In multiple regression, there was no statistical difference in the average relative importance between variables of urban, water chemistry, and flooding in explaining the observed variation in forest composition and structure. Furthermore, the combination of all three variables in the models explained 70% of the variation of the data, on average, leaving 30% to be explained by other non-measured factors. Thus, while water quality and flooding are important in the structure and composition of mangrove forests, as has been demonstrated in multiple previous studies, land cover and urban intensity may also play a role and should be considered in urban mangrove studies.

Conclusions
Urbanization presents a conservation and management challenge along tropical coastlines and especially in Caribbean mangrove systems [4,92,93]. Not only are mangrove forests lost more quickly around some of the largest cities in the world, but there is growing evidence of systemic changes to ecological function and ecosystem service provisioning in the remaining urban stands [14,43]. Using quantified metrics of urbanness, the present study corroborated these observations of novel urban mangroves by showing the most urban forests are characterized by higher maximum stem diameters, more species, and fewer truemangrove species. The inclusion of flooding and water chemistry metrics further showed that nitrogen may be an important driver of forest structure in eutrophic urban waters, that extreme high water temperatures (unrelated to urbanness) may be inhibiting stand basal area and biomass, that the relative abundance of certain species is linked to porewater salinity, which may be driven by changes in tidal connectivity along the urban gradient but that changes in hydrology due to urban infrastructure projects cannot fully account for the observed patterns in floral composition. Thus, the compounding environmental complexities of both the mangrove biome and urban landscapes create a novel suite of ecosystem drivers, the understanding of which demands well-designed and focused experiments and field studies. Furthermore, the importance of this understanding is elevated by a greater demand for valuable mangrove ecosystem services in densely populated cities. Services such as shoreline protection and fisheries support, for example, may not be well provisioned in thin urban forests dominated by Laguncularia spp. and lacking the more valuable (for these services) Rhizophora spp. and Avicennia spp. [94][95][96]. Thus, managing urban mangroves towards more sustainable social-ecological systems requires not only the conservation of remaining stands, but also a greater understanding of the influence the urban landscape has on mangrove ecology. Future work in these forests should utilize objective and quantified urban metrics for comparisons across sites, and may assess faunal communities, ecophysiology, and biogeochemistry, among others, for a more complete understanding of urban mangrove ecology and more informed management of mangroves in the Anthropocene.
was supported as part of the Next Generation Ecosystem Experiments-Tropics, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research. We are also grateful to the two anonymous reviewers who provided valuable feedback.
Funding: Work by SM was supported by an award from the US Forest Service, International Institute of Tropical Forestry (16-CA-11132762-351       Mean site forest structure and composition measurements with their strongest predictor variables from the analysis, accompanied by boxplots of least urban, urban, and most urban plots. Letters above boxes denote statistical differences. Site abbreviations used throughout the study and their corresponding locations as well select urban variables from a 1000 m sampling radius, as adapted from Branoff [45]. The urban index is a unitless metric of urbanness that incorporates several variables, with 100 being the most urban and 0 the least urban site relative to the others. Urban and mangrove cover are km 2 km −2 , street density is km km −2 , population density is people km −2 , and the ratio of urban to mangrove cover (Urban:Mang) is unitless.   Variables used to explain the variation in forest structure and composition metrics described in Table 2.   Average and standard error of stem, basal area, and biomass per hectare of all species within the three watersheds included in the study, ordered by average stem density across sites. One individual each of seven additional non-identified species were found in one plot of one site in the San Juan Bay Estuary. Halophytic plant types and salinity tolerances (ppt) provided as superscripts and obtained from Santos et al. [50], when available.  Salinity tolerance (ppt): SW seawater salinity [50].
Forests. Author manuscript; available in PMC 2021 October 21.