Effects of maternal genotypic identity and genetic diversity of the red mangrove Rhizophora mangle on associated soil bacterial communities: A field‐based experiment

Abstract Loss of plant biodiversity can result in reduced abundance and diversity of associated species with implications for ecosystem functioning. In ecosystems low in plant species diversity, such as Neotropical mangrove forests, it is thought that genetic diversity within the dominant plant species could play an important role in shaping associated communities. Here, we used a manipulative field experiment to study the effects of maternal genotypic identity and genetic diversity of the red mangrove Rhizophora mangle on the composition and richness of associated soil bacterial communities. Using terminal restriction fragment length polymorphism (T‐RFLP) community fingerprinting, we found that bacterial community composition differed among R. mangle maternal genotypes but not with genetic diversity. Bacterial taxa richness, total soil nitrogen, and total soil carbon were not significantly affected by maternal genotypic identity or genetic diversity of R. mangle. Our findings show that genotype selection in reforestation projects could influence soil bacterial community composition. Further research is needed to determine what impact these bacterial community differences might have on ecosystem processes, such as carbon and nitrogen cycling.


| INTRODUC TI ON
The world is currently experiencing losses in biodiversity at an accelerated rate (IPBES, 2019), and much research has focussed on understanding the impacts of this loss on ecosystem functioning (Duffy et al., 2017;Hooper et al., 2005;Loreau et al., 2001). As primary producers, plants play a critical role in supporting most terrestrial ecosystems, and there is a general consensus that reductions in plant diversity will lead to decreases in plant productivity (Hooper et al., 2012) with knock on effects for associated species and recycling of essential nutrients (Cardinale et al., 2011;Hooper et al., 2005). While much of the research examining plant diversity-ecosystem functioning relationships has focussed on loss of plant species or functional groups (Balvanera et al., 2006;Hooper et al., 2005;Wardle, 2016), other studies have also found that loss of genetic diversity within a dominant or foundation plant species can have similar consequences for associated species and ecosystem functioning (Crutsinger et al., 2006;Hughes et al., 2008). Genetic diversity can decline more rapidly than species diversity in response to anthropogenic pressures and habitat fragmentation, so it can be considered as a potential indicator for future species losses (Helm et al., 2009).
There can be competitive advantages for genetically diverse populations including enhanced productivity (Crutsinger et al., 2006) and reproductive success (Newman & Pilson, 2006). Genetically diverse populations have also been shown to better utilize available resources, reducing competition between conspecifics (Boyden et al., 2008), and display greater stability (Prieto et al., 2015) and capacity to adapt to environmental change or disturbance (Jump et al., 2009;Reusch et al., 2005). When genetic diversity is low, there are less likely to be individuals with the necessary traits to cope with environmental change, although there are other potential mechanisms to cope with change such as phenotypic plasticity (Gratani, 2014) and epigenetic responses (Verhoeven et al., 2016).
The level of genetic variation within a plant species can also affect interactions with other species. Crutsinger et al. (2008) found that high intraspecific genetic diversity within populations of the perennial plant Solidago altissima can deter biological invasion of other plant species. Thus, loss of genetic diversity can increase susceptibility of a population to invasion. On the other hand, genetic diversity within a parasitic species can increase their successful establishment within species-rich habitats, as has been observed in populations of the parasitic plant Rhinanthus minor (Rowntree & Craig, 2019). Changes in plant genetic diversity can also impact associated species that depend on those plants for food or shelter, with potential cascading effects throughout the food chain. For example, a meta-analysis of manipulative experiments found that richness and abundance of most arthropod trophic groups generally increase with plant genetic diversity (Koricheva & Hayes, 2018).
To date, most studies of plant genetic diversity effects on associated species have focussed on other plant species or arthropods.
There is a dearth of knowledge on how plant genetic diversity might impact associated soil microbial communities. Soil microbial communities provide many important ecosystem functions and services (Bardgett & van de Putten, 2014), and recent studies show that the structure and/or diversity of soil microbial communities significantly relates to soil functionality (Zhou et al., 2020), and processes of nitrogen (Zheng et al., 2019), and carbon cycling (Maron et al., 2018). In a grassland microcosm experiment that manipulated the soil microbiome, microbial richness and network complexity were shown to positively influence several ecosystem functions related to nutrient cycling (Wagg et al., 2019). Reduced microbial diversity is thought to lead to lower levels of soil functioning due to fewer taxa present to support the same function and lower diversity of taxa to support different functions (Wagg et al., 2019). It is therefore important to understand how changes in plant diversity may affect associated soil microbial communities and the functions that they provide.
Positive correlations have been observed between interspecific plant diversity and soil microbial richness across a range of biomes (Liu et al., 2020). Few studies have focussed on intraspecific diversity, however, in one study of Populus spp stands, relationships were detected between plant genetic diversity and soil microbial community composition, which were associated with changes in exoenzyme activity and nitrogen availability in the soil (Schweitzer et al., 2011).
In addition to genetic diversity, the identity of genotypes or phenotypes within populations of primary producers has been shown to affect ecosystem functions associated with nutrient cycling (Raffard et al., 2019). Many studies have detected differences in soil microbial community composition among plant genotypes (e.g., Aira et al., 2010;Gallart et al., 2018;Schweitzer et al., 2008;Shenton et al., 2016). Purahong et al. (2016) found effects of both genetic diversity and tree genotypic identity on soil enzyme activity in a subtropical forest, but the results were not consistent across all four tree species in the study. Plant genotypic effects on soil microbial communities have generally been attributed to differences in the chemical composition of root exudates (Micallef et al., 2009) and quantity and quality of leaf litter among genotypes (Schweitzer et al., 2005), but they could also be related to factors such as variation in root morphology or belowground carbon allocation.
In studies of natural systems, it can be difficult to separate effects of plant intra and interspecific diversity. Therefore, ecosystems that are poor in plant species could be better suited for studying intraspecific genetic diversity without the confounding effect of species diversity. Neotropical mangroves typically have few plant species, with only three true mangrove tree species found in Florida (Spalding et al., 2010), and the distribution of these vary spatially across the tidal zone (McKee, 1993). In these species-poor environments, it has been suggested that plant intraspecific genetic diversity may play an important role in shaping the communities of mangrove associate species and ecosystem functioning (Farnsworth, 1998). While we are not aware of any studies that have tested effects of intraspecific genetic diversity in mangrove ecosystems, maternal genotype has been found to influence both resistance to mortality from infestation by the parasitic beetle, Coccotrypes rhizophorae (Devlin, 2004), and seedling survival in different positions in the intertidal zone (Proffitt & Travis, 2010) in the red mangrove, Rhizophora mangle.
These results demonstrate how genotypic identity or reduction in genetic diversity could impact the ability of mangrove plant populations to deal with pressures such as herbivory or sea level rise.
Mangroves are considered to be important coastal ecosystems for carbon storage (Alongi, 2012;Donato et al., 2011), as their waterlogged soils have low rates of decomposition (Cebrian, 1999), resulting in belowground carbon stocks that are several times greater than that of other forest types (Donato et al., 2011). Human activities can contribute to the loss of plant genotypes and genetic diversity in mangroves through disturbance (Salas-Leiva et al., 2009) and fragmentation of populations (Granado et al., 2018). If loss of plant genetic diversity then affects the composition of soil microbial communities, this could have implications for soil functioning leading to changes in carbon storage rates and emissions of greenhouse gases, including carbon dioxide, nitrous oxide, and methane.
Here, our aim was to investigate how genotypic identity and genetic diversity within the mangrove tree species R. mangle affects the composition and richness of its associated soil bacterial community as determined using terminal restriction fragment length polymorphisms (T-RFLP). We hypothesized that: (a) maternal genotype would be an important determinant of soil bacterial community composition; and (b) R. mangle populations of higher genetic diversity would support richer soil bacterial communities. This was tested in a field-based experiment where plots were planted with three levels of R. mangle genotypic diversity using a pool of eight maternal genotypes. We focus on bacterial communities as they dominate the mangrove soil microbiome (Alongi, 1988;Andreote et al., 2012) and are responsible for most of the soil carbon flux in these environments (Holguin et al., 2001). We anticipate that our results will provide new insights into the influence of plant genotype and genetic diversity on belowground microbial communities in mangroves, which could help inform the restoration of degraded mangrove sites.

| Study site
The study site ( Figure 1 is sandy with shell fragments, and at the time of sampling, organic matter content in the top 10 cm of sediment in vegetated areas was around 5%-15%. Mean water levels in the Indian River Lagoon are typically higher in late autumn (Smith, 1986) when the site can be inundated for several weeks at a time, but the soil surface is generally exposed at low tide for the rest of the year.

| Experimental design
The R. mangle genetic diversity field-based experiment was set up at self-fertilizes (Lowenfeld & Klekowski, 1992) and selfing rates of 0.13-0.61 have been observed in R. mangle populations near the study area (Kennedy et al., 2017). Maternal cohort was used as a proxy for genotype; hereafter maternal cohorts are referred to as genotypes.
The experiment was set up with three levels of genetic diversity: monocultures of each of the eight genotypes, and varying mixtures of three or six genotypes, selected from a pool of the eight genotypes. The number of genotypes in the experimental pool was based on the number of trees at the collection sites that were mast producing that year (i.e., had sufficient numbers of propagules for suitable replication in the experiment). While no data are available regarding the number of genotypes that may naturally be found at the spatial scale used in this experiment, both reproduction and recruitment of R. mangle is generally low in the years following hurricanes (Proffitt et al., 2006). Thus we would not expect genotypic diversity to be as high as has been observed in similar experiments in other systems (e.g., Crutsinger et al., 2008). Experimental plots were randomly assigned to treatments and were established on unvegetated mudflats within the tidal channel. Plot placement was determined on site during planting taking into account water depth and distance from other plots (at least 2 m). Each plot measured 45 x 30 cm and 12 propagules were randomly assigned within a matrix of three rows with 15 cm spacing between each propagule (based on an average distance observed between propagules, which had naturally recruited at a nearby site). Propagules that failed to establish with root growth 1 week after planting were replaced by another propagule of the same cohort. Some maternal trees yielded more propagules than others, so genotypes were assigned to plots in a stratified random design taking number of available propagules per genotype into account. This resulted in uneven representation of the different genotypes. The number of replicate monoculture plots per genotype ranged from 4 to 6. In the three genotype treatment, each genotype was represented in 9-17 plots. In the six genotype treatment, each genotype was represented in 6-13 plots. There were 81 plots in total: 39 plots were planted as monocultures (12 propagules from one genotype), 29 plots were planted with three genotypes (four propagules from each), and 13 plots were planted with six genotypes (two propagules from each). With a pool of eight genotypes, there are more possible combinations of three genotypes compared to six genotypes, hence why more three genotype plots were planted. In June 2016, 1 month prior to our soil sampling, mean tree height within experimental plots was 1.9 m and 88% of trees initially planted had survived. All trees within experimental plots had produced multiple prop roots and most trees had reached reproductive maturity.

| Soil sampling
In August 2016, almost five years after the experiment was set up, we sampled soil from 47 of the original experimental plots (Figure 1). The samples were collected in the summer as this is when we would expect the bacterial community to be most active at this site due to warm temperatures and increased exposure to oxygen (compared to late autumn/winter when the soil surface is not exposed at low tide). To test for genotype effects on soil bacterial communities, we sampled 24 of the monoculture plots (three replicates of each genotype). To test for genotypic diversity effects, we also sampled 12 plots that had been planted with three genotypes, and 11 plots that had been planted with six genotypes. We only sampled from plots whose identity could be confirmed by above ground tags in the field. Some tags had become buried in the time since planting and excavation of the below ground tags would have caused significant disturbance to the soil and associated microbial communities. Six subsample soil cores were collected in each plot from random locations among the R. mangle prop roots, to a depth of approximately 10 cm, using corers made from sterile 50 ml syringes with the end cut off. The cores for each plot were placed in a single ziplock bag and mixed to form a composite sample of bulk soil. Five control soil samples (also consisting of six subsamples) were collected from mudflat areas across the site that were at least 2 m from the nearest visible tree prop root or pneumatophore. Soil layers below 10 cm were not sampled as oxygen concentrations and hence, microbial activity, rapidly decline with depth from the surface in mangrove soils (Booth et al., 2019). Pore water samples were extracted under suction from each plot (and control location) using a sipper (McKee et al., 1988) and collected into 25 ml scintiallation vials. The GPS coordinates and number of remaining trees were recorded for each sampled plot. All samples were collected during morning low tide and were kept on ice in a cool box in the field and during transportation to the laboratory.

| Soil and bacterial analyses
Pore water pH and salinity were measured on the day of collection, after the pore water had equilibrated to room temperature, using an ExStik pH meter (Extech) and an ORA 1SB analog refractometer (Kern Optics) respectively. A subsample from each composite soil sample was oven-dried at 60°C for 48 hr then ground in a Mixer Mill MM400 (Retsch GmbH). Total % carbon and nitrogen were measured by dry combustion of 10 mg of the ground soil using a Vario Negative PCR controls were also prepared in the same way except instead of adding DNA from the soil samples, 1 µl of sterile water was added. PCR amplification was carried out in a SimpliAmp thermal cycler (Thermo Fisher Scientific) under the following conditions: 94°C for 1 min 30 s; followed by 30 cycles of 94°C for 45 s, 55°C for 1 min and 72°C for 3 min; followed by 10 min extension at 72°C. PCR products were purified using ethanol precipitation and resuspended in 15 µl molecular biology grade water. 5 µl purified PCR product was digested for 2 hr at 37°C using the MspI enzyme (Promega). Samples were purified again using ethanol precipitation and resuspended in 10 µl molecular biology grade water. 1 µl of the final sample was then mixed with 8.6 µl Hi-Di formamide and 0.4 µl GeneScan 600 LIZ size standard for fragment sizing using an ABI 3730 (Applied Biosystems).
The resulting data were analyzed using GeneMapper version 4.1 (Applied Biosystems) to score peaks over 150 fluorescence units in the range 50-500 bp and assign bins. This process resulted in a binary matrix file with presence/absence of each bin within each sample. The peak profile of each sample was visually checked and manually edited where peaks were observed that were not detected by the software.
All DNA samples were run in duplicate PCRs and only bins that appeared in both replicates were included in the final binary matrix. Only bins that occurred in two or more samples were included in analyses.
To compare soil bacterial community profiles, Jaccard similarity matrices were constructed and visualized with nonmetric multidimensional scaling (NMDS) using the metaMDS function. The adonis function was used to perform permutational multivariate analysis of variance (PERMANOVA) with 1,000 permutations to test for significant effects of genotypic identity or genetic diversity on bacterial community composition. As terms are added sequentially, all measured environmental variables (total soil nitrogen and carbon, porewater pH and salinity, and number of trees remaining in the plot) were included in the model before the treatment so that variability attributed to the environmental variables was removed before testing for treatment (genotypic identity or genetic diversity) effects. Environmental vectors were fitted to the NMDS using the envfit function. When testing for genotypic identity, only monoculture plots were included in the analysis and genotypic identity was nested within origin site of the propagules. Mudflat soil samples were not included in the PERMANOVAs testing for genotypic identity or genetic diversity effects. Homogeneity of group dispersion (variance) was measured using the betadisper function and differences between groups were tested using the permutest.betadisper function.
The number of T-RFLP bins per sample was used as a proxy for bacterial taxa richness. Kruskal-Wallis tests were used to test whether bacterial taxa richness or any of the measured environmental variables differed between the treatments. When significant differences were found, post hoc Dunn tests were performed using the dunnTest function in the FSA package (Ogle et al., 2020) with p-values adjusted for multiple testing using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995).

| RE SULTS
Due to the hydrological flow at the site, gradients of pore water salinity, total soil nitrogen, and total soil carbon were observed across plots in the interpolation maps ( Figure S1). However, the random distribution of the experimental plots meant that these variables did not significantly differ among the diversity treatments or among the genotypes in the monoculture plots when tested with Kruskal-Wallis tests (Table 1); as such, these gradients should not confound our results. Within the treatment plots, pore water salinity ranged from 35‰ to 55‰, total soil carbon ranged from 3.97% to 8.84%, and total soil nitrogen ranged from 0.07% to 0.49%. There was no obvious pattern of pore water pH across the site, which ranged from 6.70 to 7.59.
An average of 10 trees remained across the sampled treatment plots (range 5-12 trees) and there was no difference in propagule survival between the treatments at this stage of the experiment (Table 1). Due to the amount of plant growth since the experiment had been setup and the difficulty of locating ID tags, we were unable to identify which of the original propagules were no longer present.
Thus we could not confirm the diversity level remaining in all plots at the time of sampling. None of the environmental variables measured significantly differed between the mudflat and the treatment plots.

| Genotypic identity
Within the monoculture plots, PERMANOVA detected a significant relationship between the maternal genotypic identity of R. mangle (when nested within origin site) and the associated soil bacterial community (Table 2), with genotypic identity explaining 34% of the variation in bacterial community composition. This was independent of origin site of the propagules, which was not found to be a significant factor explaining bacterial community composition. The ordination plot (Figure 2) reveals that while some of the genotypes cluster quite closely together, indicating that they have similar soil bacteria communities, a few genotypes occupied space that was distinct from some of the other genotypes. Bacterial community composition within the monoculture plots was also significantly related to total soil nitrogen ( Table 2). Bacterial taxa richness and the measured environmental variables did not differ among genotypes (Table 1), and there was no detectable difference in group dispersion of soil bacterial communities between genotypes (F 7,16 = 0.957, p = 0.535).

| Genetic diversity
Tree presence had a clear effect on soil bacterial community composition with mudflat samples clustering apart from the samples collected within the treatment plots ( Figure 3). However, within the experimental plots, there was no significant difference among the R. mangle genetic diversity treatments (Table 3). Of the environmental variables measured, only total soil nitrogen was significantly related to bacterial community composition (Table 3). Although the six genotype treatment appeared to be more variable, there was no significant difference in group dispersion between the genetic diversity treatments (F 2,44 = 2.065, p = 0.169; Figure 2).
Bacterial taxa richness trended toward a humped distribution, increasing from monocultures to the three genotype treatment, then falling back in the six genotype treatment (Figure 4), but this relationship was not significant (Table 1). The measured environmental variables did not differ between genetic diversity treatments (Table 1). When the mudflat samples were included in the analysis, there was an effect of treatment group (χ 2 = 11.66, df = 3, p = 0.009) with taxa richness being lower in the mudflat samples compared to the one (p = 0.016) and (p = 0.005) three genotype treatments.

| D ISCUSS I ON
Our study set out to determine whether genotypic identity and genetic diversity of the neotropical red mangrove, R. mangle, influence associated soil bacterial communities. We demonstrate that maternal genotypic identity can drive changes in associated soil bacterial community composition in wetland environments. If different genotypes host different soil bacterial communities, it is plausible that plant community assemblages consisting of more genotypes (i.e., that are more genetically diverse) would harbor more diverse bacterial communities. However, this was not observed in our study, in that richness of the bacterial communities did not differ between the R. mangle genetic diversity treatments.

Our findings support those of a growing number of studies
showing that plant genotype affects soil microbial community composition (e.g., Aira et al., 2010;Gallart et al., 2018;Schweitzer et al., 2008;Shenton et al., 2016

TA B L E 2 Permutational multivariate analysis of variance (PERMANOVA)
test results for effects of Rhizophora mangle maternal genotypic identity, maternal origin site, and environmental variables within monoculture plots on the Jaccard similarity index of soil bacterial communities based on 16S T-RFLP profiling and growth form (Zogg et al., 2018), as well as among S. alterniflora plants from different populations grown in a common garden (Nie et al., 2010). The results from our study show that even in environments where tidal waters could promote mixing of microbial communities across a site, tree genotype still exerts a detectable influence on the composition of bacterial communities in bulk soils.
In wetlands, the rhizosphere could be particularly important in controlling bacterial community diversity and composition. Higher bacterial diversity observed in S. alterniflora rhizosphere compared to bulk soil has been attributed to the plant roots providing oxygen and labile organic compounds in an environment lacking both (Zogg et al., 2018). Our finding that bacterial taxa richness was higher in experimental plots compared to the surrounding mudflats lends support to this suggestion. The root morphology of R. mangle is highly adapted for survival in anoxic soil environments and their prop roots can transport oxygen belowground to the rest of the root system (Scholander et al., 1955). In such a plant, any differences in root morphology among genotypes could be expected to have a substantial influence on the composition of associated soil bacterial communities by controlling the supply of oxygen and exudates to the rhizosphere. Variations among genotypes in the chemical composition of exudates (Micallef et al., 2009), or the quantity and quality of leaf litter (Schweitzer et al., 2005), could also drive changes in the bacterial community.
In contrast to Schweitzer et al.'s (2011) study, we did not find an effect of plant genetic diversity on associated bacterial community composition. Schweitzer et al. (2011) found that natural stand genetic diversity was linked to microbial community composition as well as measures related to microbial soil functioning including exoenzyme activity and inorganic nitrogen concentrations. Nitrogen cycling has also been shown to be affected by plant genetic diversity in a pot experiment through increased root litter nitrogen content and decomposition rates compared with litter produced by sibling groups (Semchenko et al., 2017). Although our study and Schweitzer et al.'s (2011) did not find any relationships between plant genetic diversity and total soil nitrogen content, Schweitzer et al. (2011) detected unimodal relationships between genetic diversity and soil nitrate and ammonium concentrations (which we did not measure).
While the trees in the experiment plots had reached maturity by the time of our soil sampling, it is possible that a longer time scale is needed before small effects on nutrient cycling become detectable.
We suggest that further investigations are warranted to evaluate the relationships between plant genetic diversity, soil microbial community composition, and processes of nitrogen cycling across different biomes in field conditions.
We did not observe any differences in bacterial taxa richness between R. mangle genotypes or genetic diversity treatments. While a more diverse plant community might be expected to support a more diverse soil bacterial community, evidence for this is lacking, suggesting that other environmental factors may have a stronger role in regulating soil bacterial diversity (Fierer & Jackson, 2006;Prober et al., 2015;Wardle, 2006). However, the presence of R. mangle was found to have some influence on bacterial taxa richness, with significantly higher richness found in the treatment plots compared to the surrounding mudflats. Our bacterial community profiles were based on presence/absence data, so we were only able to assess richness of T-RFLP fragments as a measure of diversity. Next generation sequencing could be used to determine abundance and evenness of different bacterial taxa in order to calculate diversity indices, which could provide a better understanding of potential effects of plant genetic diversity on bacterial community diversity. Measures of soil function, such as rates of respiration, mineralization, and decomposition, or use of qPCR or metagenomics to assess abundance changes in functional genes could also be used to ascertain whether nutrient F I G U R E 2 NMDS ordination plot based on the Jaccard's similarity index of bacterial communities examined by 16S T-RFLP in soil samples collected within monoculture plots (stress = 0.121). Each color represents a different Rhizophora mangle maternal genotype and each triangle corner is one plot of that genotype F I G U R E 3 NMDS ordination plot based on the Jaccard's similarity index of bacterial communities examined by 16S T-RFLP in soil samples collected from plots of Rhizophora mangle with different levels of maternal genotypic diversity and control samples collected from the surrounding mudflats (stress = 0.137). The only fitted environmental vector found to be significantly related to soil microbial community composition (total soil N) is shown with the gray arrow cycling processes are impacted as a result of changes in bacterial community composition.
Our experiment comprised maternal cohorts as a proxy for genotype, and we did not conduct genotype analyses on the propagules or the maternal trees. As such, it is likely that there is a degree of variation within each maternal cohort in addition to variation among cohorts. However, due to self-fertilization rates of R. mangle (Kennedy et al., 2017) in this region, we would expect variation to be greater among maternal cohorts than within. As we do not know how related each of our maternal cohorts are, some combinations of genotypes within the three genotype and six genotype diversity treatments may be less diverse than others. This may have impacted our analyses on effects of genotypic diversity on the soil bacterial community composition. Nevertheless, and despite the limitations of our methodology and frequent mixing of water across the site due to tidal influences, we were still able to detect differences in bacterial communities of bulk soil between R. mangle genotypes, which points to the strong effect plants can have on the composition of associated bacterial communities.
As the importance of mangroves is becoming more widely recognized, there has been an increase in mangrove restoration efforts in many regions (Bosire et al., 2008 (Hines et al., 2014), further work is needed to understand the potential implications of genetic diversity loss on important soil functions associated with carbon and nutrient cycling. Changes to the rate of decomposition of organic matter in these environments could affect not only their carbon storage capacity, but also their ability to keep up with sea level rise through vertical accretion of the soil.

ACK N OWLED G M ENTS
The study contained in this publication was supported by the Our thanks also go to three anonymous reviewers for taking the time to provide helpful comments to improve the quality of this publication.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data associated with this manuscript is available from the Dryad Digital Repository (https://doi.org/10.5061/dryad.5tb2r bp2x).