Rhizophora is the dominant genus of mangrove forests on the Atlantic coast of northern South America. What determines the zonation frequently observed in sympatric populations of the two neotropical species R. mangle, R. racemosa, and their hybrids, R. × harrisonii, is an open question. The most widely held hypothesis is that differences in salinity tolerance among the taxonomic groups explains the observed zonation. To address this question, we analyzed the elemental composition of soils and canopy leaves from 60 Rhizophora spp. trees distributed in different intertidal zones of an estuarine site of the Paria Gulf, Venezuela. The low intertidal zone showed lower salinity, organic matter, C, N, S, and Na, and higher bulk density, Al, Fe, and Mn concentrations compared with the higher intertidal zones. Using morphological characters and microsatellite molecular markers, we identified 39 pure R. mangle, 19 hybrids, and only two pure R. racemosa. We found that both intertidal position and taxonomic groups explained most of the differences in leaf variables measured across trees. The ratio Mg/Ca, however, was higher in R. mangle than in R. racemosa and hybrids regardless of intertidal position. Moreover, at some specific intertidal position, R. mangle differed from R. racemosa and hybrids in the values of C, N, K, Mg, Fe, Mn, C/N, K/Ca, S/Ca, and d13C. We conclude that despite the scarcity of R. racemosa and the absence of a clear species zonation, our results suggest that R. mangle copes with salinity differently than R. racemosa and R. × harrisonii.

The study site is located in an estuarine-riverine mangrove forest dominated by red (Rhizophora spp.) and black (Avicennia germinans) mangroves with few individuals of white mangrove (Laguncularia racemosa). The location of this mangrove forest is within the Paria gulf in a site known as Punta de Piedra (10.421767 N; 62.796708 W), bordering Turuépano National Park (Colonnello et al., 2009). Tidal channels such as Turuépano, La Palma, Guariquén, Guarapichito, and Ajíes surround this area. In addition, the San Juan River discharge is 25 km southeast of the site (Lasso et al., 2004), so freshwater and sediment discharge strongly influence the mangrove forest of the area. 

The study site was selected based on its accessibility, density, and homogeneity of the vegetation, and relative protection from the Paria Gulf heavy winds. Based on information from the closest meteorological stations (i.e., Caripito and Güiria), the area has a short dry season between February and March (Walter & Medina, 1971). However, data from a climate model covering a more extensive period (1982-2012) indicate a longer dry season from January to April, with averages of 1444 mm of precipitation and 26.9ºC for temperature (Climate-Data.org, 2016). The river freshwater discharge and the effect of large tidal amplitudes (≈ 1 m, Herrera et al., 1981) compensate for the reduction in rainfall during the first quarter of the year. From June to December, rainfall exceeds 100 mm per month, contributing to significant leaching of salt from superficial soil layers. River discharge and rainfall maintains water salinity levels well below standard seawater (35‰) at the fringe of the mangrove-lined coast. For this study, we conducted fieldwork in December 2010 at the beginning of the dry season. 

Delimitation of the intertidal zones and sampling design

The tidal regime in the Paria Gulf is bi-diurnal and the tidal range varies between 1.8 and 3.2 m. Therefore, the penetration of tidal waters into the forest varies accordingly. In this study area, the maximum penetration inland of the mangrove forest is about 200 m. Beyond 200 m, non-halophytic species such as Symphonia globulifera and Cassipourea spp. (Colonnello et al. 2009) begin to appear.

Based on this characterization of the study site, we established two transects perpendicular to the coast border toward the upper tidal zone and into the mangrove forest, separated 200 m from each other. On each transect, we delimited the two extremes of the intertidal range. The low intertidal zone (LI) was delimited as a 20m-wide stripe located at 2-20 m from the shoreline. Similarly, the upper intertidal zone (UI), was delimited as a 20m-wide strip at 150-170 m from the shoreline. Finally, a 20m-wide strip at 90-110 m from the shoreline was considered to represent the middle intertidal zone (MI). This intermediate area helped us test for gradients of different parameters across the intertidal range. Thus, our sampling design included two plots in each of the three intertidal zones, one in each transect, totaling six plots.

At the low intertidal zone (LI), each one of the two plots had an area of 900 m2 (30 x 30 m) totaling 1,800 m2. Because of the density of Rhizophora spp. trees with DBH > 5 cm increased from the shoreline (LI) toward the middle (MI) and upper (UI) intertidal zone, we reduced the size of the plots progressively. At the MI zone, each one of the two plots covered 400 m2 (20 x 20 m), totaling 800 m2. At the UI zone, each one of the two plots covered 225 m2 (15 x 15 m), totaling 450 m2. 

Measurements of interstitial water salinity

We measured interstitial water salinity in situ at eight distances, along the two transects, from the coastal border toward the upper tidal zone into the mangrove forest, spanning the three intertidal zones. The distances were 2 and 20 m (at LI); 90, 100, and 110 m (at MI); and 150, 160, and 170 m (at UI). At each distance point, we collected porewater samples at 20, 40, and 80 cm depth with a 1 m long aluminum tube with perforations at the end, connected to a 60-mL syringe. Salinity was measured with a digital refractometer (PAL-06S, Atago) and expressed in ‰. We used one or two replicates for each depth for a total of 72 measurements across the study site.  

​​​​​​​Sampling of soils and canopy leaves of Rhizophora spp. individuals

We randomly collected four soil samples at 0-20 cm depth of each plot for 24 samples (i.e., four random samples x three intertidal zones x two plots per intertidal zone). We discarded one sample during processing. The remaining 23 samples were dried in the laboratory using a ventilated oven at 45oC to constant weight, ground, and then passed through a 20-mesh sieve. We separated each sample into two subsamples to measure the elemental composition (see section 2.5) and pH. For pH, the subsamples were hydrated with deionized water (1:5), homogenized manually for 5 minutes, and allowed to rest for 14 h before measurements with a HANNA instrument (model 991301). Following the same sampling design, we measured bulk density in soil samples collected between 0-10 cm depth using PVC tubing (5 cm diameter, 10 cm length) for 24 samples. Samples were dried in the laboratory using a ventilated oven (45oC) until constant weight. This measurement was the quotient of soil dry mass/original volume of the PVC cylinder (200 cm3).

Within each one of the six plots, we randomly selected ten trees of Rhizophora spp. with inflorescences totaling 60 individuals. Eight fully expanded green canopy leaves were collected from each of these trees for laboratory analysis of elements and isotopes (see sections 2.5 and 2.6).

​​​​​​​Elemental composition of canopy leaves and soils

Sixty samples of canopy leaves were oven-dried at 65oC, grounded and subsequently passed through an 18-mesh sieve. The powdered samples were ashed at 490oC and then dissolved in HCl 6N.  Using the thermogravimeter (LECO TGA701), we measured dry mass (105oC) and ash content (490o).  In parallel, the 23 soil samples were digested sequentially in concentrated hot HNO3, H2O2, and HCl (Huang & Schulte, 1985). Then, using digested subsamples of both leaves and soil samples, Na, K, Mg, Ca, Fe, Mn, S, and P concentrations were measured using the Spectro Analytical Instruments GmbH (Boschstr. 10, 47533 Kleve, Germany).  Additionally, total carbon and nitrogen were also measured using a macro dry combustion method with the Leco CNS- 2000 Analyzer (Leco Corp. 3000 Lakeview Ave. St. Joseph, Michigan 49085, USA).

​​​​​​​Leaf isotopic C and N ratios

Isotopic C and N ratios of the leaves sampled were measured for each one of the 60 Rhizophora spp. trees at the L7 isotope laboratory of the University of Miami, Coral Gables (http://penguin.bio.miami.edu). Small tin spheres containing samples were placed in an automated elemental analyzer (Eurovector, Milan, Italy) and pyrolyzed. Gases from the pyrolysis were led into a mass spectrometer (Isoprime, Elementar) and analyzed for 13C and 15N abundance. Carbon and nitrogen isotope ratios are expressed as  d13C or  d15N = [(Rsample/Rstandard) -1] ·103. Rsample and Rstandard represent the molar ratio of the heavier (13C and 15N) to lighter isotope (12C and 14N) in the sample and standard, respectively. The standard for carbon isotope ratio was belemnite from the PeeDee formation in South Carolina, with an analytical precision of 0.1 ‰. The standard for the nitrogen isotope ratio was atmospheric nitrogen, with an analytical precision of 0.3 ‰.

​​​​​​​Mangrove species identification by morphological characteristics

In each of the six plots at three different intertidal zones, we identified mangrove trees of Rhizophora spp., Avicennia germinans, and Laguncularia racemosa. Identification of A. germinans and L. racemosa was straightforward because there are no other species of these genera in the area; however, identification of Rhizophora trees was complex. The only morphological character distinguishing the two species and the hybrids is the inflorescence type (i.e., the number of bifurcations in the inflorescence). Based on the inflorescence type, the 60 sampled Rhizophora trees were identified as follows: those with 1-2 bifurcations as R. mangle, those with 5-6-7 bifurcations as R. racemosa, and those with 3-4-5 bifurcations as hybrids, R. × harrisonii (Jiménez 1987; Cerón-Souza et al., 2010). Due to the partly overlapping ranges of variation of inflorescence bifurcation in the three taxonomic units, which is associated to ongoing hybridization and introgression, field identification should ideally be corroborated genetically with molecular markers (Cerón-Souza et al., 2010; Cerón-Souza et al., 2014). 

​​​​​​​Genetic identification of sampled Rhizophora individuals

Of the 60 Rhizophora trees sampled, we genotyped 46 individuals to validate the taxonomic identification based on inflorescence type using DNA extracted from leaves. We used ten microsatellite loci (i.e., RM07, RM11, RM19, RM21, RM36, RM46, RM05, RM50, RM59, RS67) following the protocols in Cerón-Souza et al., (2010). The genetic composition of individuals and the hybrid class were determined using Bayesian genetic clustering analyses with Structure 2.3.3 (Falush et al., 2003, 2007; Pritchard et al., 2000) and NewHybrids 1.1 beta (Anderson & Thompson, 2002). In neither of the two analyses was the taxonomic identification based on inflorescence-types used as prior information. For Structure, an admixture model was assumed, with correlated allele frequencies and a uniform prior probability of K. We ran ten iterations for each number of groups assayed, from K = 1 through K = 10, using 500,000 replicates of MCMC after a burn-in of length 50,000 replicates. We applied the criterion of Evanno et al. (2005) to select the best estimate of the number of genetic groups based on DK calculation. For NewHybrids, we performed the analysis five times, each starting with a different random number of seeds and including 500,000 iterations and 50,000 burn-in steps of the MCMC chain using a uniform allele and genetic frequency prior. Based on this approach, the 46 individuals were assigned to one of six possible genotype classes: 1) pure R. mangle parent, 2) pure R. racemosa parent, 3) F1 hybrid (50 % of the genome originated from R. mangle and 50 % from R. racemosa), 4) F2 hybrid (50 % from F1 hybrid and 25 % from each of the parents R. mangle and R. racemosa), 5) backcrosses with R. mangle (50 % from F1 hybrid and 50 % from parent R. mangle), and 6) backcrosses with R. racemosa (50 % from F1 hybrid and 50 % from parent R. racemosa).

​​​​​​​​​​​​​​Statistical analysis

To test differences in soil characteristics, we used a one-way ANOVA comparing LI, MI, HI intertidal zones. We also performed a principal component analysis of element concentration and element ratios to determine their correlation and characterize their spatial distribution across the intertidal zone. For the interstitial soil salinity, we used a two-way ANOVA to compare the effect of intertidal zone and depth.

To test element concentrations, ratios of elements, and isotope differences across Rhizophora individuals, we performed a two-way ANOVA comparing the effect of the position in the intertidal zone and the taxonomic identification combining morphology and molecular markers. 

In all analyses for both soils and leaves variables, we confirmed the homogeneity of variances (Bartlett test). For the few variables in which unequal variances were detected, we used Welch's test for means allowing for unequal standard deviations. We compared the means using the Tukey-Kramer HSD test at p<0.05. For all statistical tests, we used JMP 13 (SAS Institute Inc. 2016).