In a time of global change, having an understanding of the nature of biotic and abiotic factors that drive a species’ range may be the sharpest tool in the arsenal of conservation and management of threatened species. However, such information is lacking for most tropical and epiphytic species due to the complexity of life history, the roles of stochastic events, and the diversity of habitat across the span of a distribution. In this study, we conducted repeated censuses across the core and peripheral range of Trichocentrum undulatum, a threatened orchid that is found throughout the island of Cuba (species core range) and southern Florida (the northern peripheral range). We used demographic matrix modeling as well as stochastic simulations to investigate the impacts of herbivory, hurricanes, and logging (in Cuba) on projected population growth rates (? and ?_s) among sites.

Field methods

Censuses took place between 2013 and 2021. The longest census period was that of the Peripheral population with a total of nine years (2013–2021). All four populations in Cuba used in demographic modeling that were censused more than once: Core 1 site (2016–2019, four years), Core 2 site (2018–2019, two years), Core 3 (2016 and 2018 two years), and Core 4 (2018–2019, two years) (Appendix S1: Table S1). In November 2017, Hurricane Irma hit parts of Cuba and southern Florida, impacting the Peripheral population.

The Core 5 population (censused on 2016 and 2018) was small (N=17) with low survival on the second census due to logging. Three additional populations in Cuba were visited only once, Core 6, Core 7, and Core 8 (Table 1). Sites with one census or with a small sample size (Core 5) were not included in the life history and matrix model analyses of this paper due to the lack of population transition information, but they were included in the analysis on the correlation between herbivory and fruit rate, as well as the use of mortality observations from logging for modeling. All Cuban sites were located between Western and Central Cuba, spanning four provinces: Mayabeque (Core 1), Pinar del Rio (Core 2 and Core 6), Matanzas (Core 3 and Core 5), and Sancti Spiritus (Core 4, Core 7, Core 8).

At each population of T. undulatum presented in this study, individuals were studied within ~1-km strips where T. undulatum occurrence was deemed representative of the site, mostly occurring along informal forest trails. Once an individual of T. undulatum was located, all trees within a 5-m radius were searched for additional individuals. Since tagging was not permitted, we used a combination of information to track individual plants for the repeated censuses. These include the host species, height of the orchid, DBH of the host tree, and hand-drawn maps. Individual plants were also marked by GPS at the Everglades Peripheral site. If a host tree was found bearing more than one T. undulatum, then we systematically recorded the orchids in order from the lowest to highest as well as used the previous years’ observations in future censuses for individualized notes and size records. We recorded plant size and reproductive variables during each census including: the number of leaves, length of the longest leaf (cm), number of inflorescence stalks, number of flowers, and the number of mature fruits. We also noted any presence of herbivory, such as signs of being bored by M. miamensis, and whether an inflorescence was partially or completely affected by the fly, and whether there was other herbivory, such as D. boisduvalii on leaves.

We used logistic regression analysis to examine the effects of year (at the Peripheral site) and sites (all sites) on the presence or absence of inflorescence herbivory at all the sites. Cross tabulation and chi-square analysis were done to examine the associations between whether a plant was able to fruit and the presence of floral herbivory by M. miamensis. The herbivory was scored as either complete or partial.

During the orchid population scouting expeditions, we came across a small population in the Matanzas province (Core 5, within 10 km of the Core 3 site) and recorded the demographic information. Although the sampled population was small (N = 17), we were able to observe logging impacts at the site and recorded logging-associated mortality on the subsequent return to the site.

Matrix modeling

Definition of size-stage classes

To assess the life stage transitions and population structures for each plant for each population’s census period we first defined the stage classes for the species. The categorization for each plant’s stage class depended on both its size and reproductive capabilities, a method deemed appropriate for plants (Lefkovitch 1965, Cochran and Ellner 1992). A size index score was calculated for each plant by taking the total number of observed leaves and adding the length of the longest leaf, an indication of accumulated biomass (Borrero et al. 2016). The smallest plant size that attempted to produce an inflorescence is considered the minimum size for an adult plant. Plants were classified by stage based on their size index and flowering capacity as the following: (1) seedlings (or new recruits), i.e., new and small plants with a size index score of less than 6, (2) juveniles, i.e., plants with a size index score of less than 15 with no observed history of flowering, (3) adults, plants with size index scores of 15 or greater. Adult plants of this size or larger are capable of flowering but may not produce an inflorescence in a given year. The orchid’s population matrix models were constructed based on these stages. 

In general, orchid seedlings are notoriously difficult to observe and easily overlooked in the field due to the small size of protocorms. A newly found juvenile on a subsequent site visit (not the first year) may therefore be considered having previously been a seedling in the preceding year. In this study, we use the discovered “seedlings” as indicatory of recruitment for the populations. Adult plants are able to shrink or transition into the smaller juvenile stage class, but a juvenile cannot shrink to the seedling stage.

Matrix elements and population vital rates calculations

Annual transition probabilities for every stage class were calculated. A total of 16 site- and year-specific matrices were constructed. When seedling or juvenile sample sizes were < 9, the transitions were estimated using the nearest year or site matrix elements as a proxy. Due to the length of the study and variety of vegetation types with a generally large population size at each site, transition substitutions were made with the average stage transition from all years at the site as priors. If the sample size of the averaged stage was still too small, the averaged transition from a different population located at the same vegetation type was used. We avoided using transition values from populations found in different vegetation types to conserve potential environmental differences. A total of 20% (27/135) of the matrix elements were estimated in this fashion, the majority being seedling stage transitions (19/27) and noted in the Appendices alongside population size (Appendix S1: Table S1). The fertility element transitions from reproductive adults to seedlings were calculated as the number of seedlings produced (and that survived to the census) per adult plant.

Deterministic modeling analysis

We used integral projection models (IPM) to project the long-term population growth rates for each time period and population. The finite population growth rate (?), stochastic long-term growth rate (?_s), and the elasticity were projected for each matrices using R Popbio Package 2.4.4 (Stubben and Milligan 2007, Caswell 2001). The elasticity matrices were summarized by placing each element into one of three categories: fecundity (transition from reproductive adults to seedling stage), growth (all transitions to new and more advanced stage, excluding the fecundity), and stasis (plants that transitioned into the same or a less advanced stage on subsequent census) (Liu et al. 2005). Life table response experiments (LTREs) were conducted to identify the stage transitions that had the greatest effects on observed differences in population growth between select sites and years (i.e., pre-post hurricane impact and site comparisons of same vegetation type).

Due to the frequent disturbances that epiphytes in general experience as well as our species’ distribution in hurricane-prone areas, we ran transient dynamic models that assume that the populations censused were not at stable stage distributions (Stott et al. 2011). We calculated three indices for short-term transient dynamics to capture the variation during a 15-year transition period: reactivity, maximum amplification, and amplified inertia. Reactivity measures a population’s growth in a single measured timestep relative to the stable-stage growth, during the simulated transition period. Maximum amplification and amplified inertia are the maximum of future population density and the maximum long-term population density, respectively, relative to a stable-stage population that began at the same initial density (Stott et al. 2011). For these analyses, we used a mean matrix for Core 1, Core 2 Core 3, and Core 4 sites and the population structure of their last census. For the Peripheral site, we averaged the last three matrices post-hurricane disturbance and used the most-recent population structure. We standardized the indices across sites with the assumption of initial population density equal to 1 (Stott et al. 2011).  Analysis was done using R Popdemo version 1.3-0 (Stott et al. 2012b).

Stochastic simulation

We created matrices to simulate the effects of episodic recruitment, hurricane impacts, herbivory, and logging (Appendix S1: Table S2). The Peripheral population is the longest-running site with nine years of censuses (eight transitions) which we used to select matrix elements that contained the years that experienced episodic recruitments, direct impact from Hurricane Irma (category 3, a major hurricane), as well as leaf herbivory impacts from scale insects. Trichocentrum undulatum experiences infrequent recruitment, but we captured episodes of high recruitment at the Peripheral site during the first and sixth monitoring years.

Specifically, an episodic recruitment simulation matrix for one of the Cuban field sites with the longest-running census period (Core 1 population) was created to calculate the probability of episodic recruitment that the population will need in order to reach a stable ?_s = 1 (Appendix S1: Table S2). The simulated matrix for Core 1 was created using the average transitions for the real-time site censuses and substituting the transition elements of fecundity, seedling survival, and seedling growth transition from the episodic recruitment event at the Everglades.

Due to D. boisduvalii scale being present during all censuses in the Peripheral population, with the exception of one year, the best way to simulate the effects of herbivory is to remove the associated mortality. To simulate the impact of the lethal D. boisduvalii attacks, we removed the plants that were observed as dead due to scale infestation during the third and fourth census and constructed a new matrix to calculate ?s (Appendix S1: Table S2). Similarly, to infer the impact of Hurricane Irma, we excluded mortality caused by the hurricane (2017) at the Peripheral site, i.e., deaths attributed to the hurricane for each stage class documented in the fifth census and generated a non-hurricane affected matrix to calculate ?_s.

Our censuses for the Peripheral population began in 2013, eight years after the last hurricane impacts in 2005 and four years before Hurricane Irma. To project the effects of hurricanes at the Peripheral site, positions of hurricanes passing through a 100 km circle centered at the site between the years 1970 and 2020 were obtained from National Oceanic and Atmospheric Administration (NOAA) and the National Hurricane Center’s HURDAT2 climatology (Landsea and Franklin 2013). We used the historical hurricane frequency at the Peripheral site (one hurricane per decade) to calculate the ?s under the same hurricane regime. In an abundance of caution for the potential increase in TC intensification into hurricanes at the site, we also carried out a simulation using a scenario with more frequent hurricane episodes (an increase in hurricane frequency of 50%) to see the effect on the projected long-term population growth of the population (Appendix S1: Table S1). We created a Markovian chain (adapted code from S. Elsner 2008) with four major transition phases gathered from the Peripheral population censuses: (1) phase I, the hurricane year (census 5); (2) phase II, first year post-hurricane (census 6); (3) phase III, second and third year post hurricane (7 and 8); (4) phase IV, non-hurricane affected years (census 1, 2, 3, and 4) (Appendix S1: Table S1). To simulate the historic and increased probability of hurricane impacts on long-term population growth rates (λ values) at the Everglades National Park we used two Markovian chains. The probability of a hurricane happening on any given year was 0.1. If a hurricane did occur, then Phase II was followed by Phase III, unless another hurricane. The two matrices in phase III occurred at equal probability on the second and third year post-hurricane. On the fourth year after a hurricane and until the next hurricane occurs, the four matrices in phase IV occurred at equal probability. To project the effects of an increase in hurricane frequency, we applied changes to the yearly hurricane probability to 0.15 (an increase of 50%). The remaining probabilities for the above post-hurricane stayed the same.

The Core 3 site is one of the largest populations censused in this study and is found in a coastal forest that experiences periodic flooding similarly to the Peripheral site. Hurricane Irma swept over the northern coast of Cuba and did not affect Core 3. We used the mortality for each stage class observed at the Peripheral population from hurricane damage and created a hurricane matrix for the Core 3 population. We simulated different hurricane year probabilities alongside the other empirical-data-based matrices for the Core 3 censuses until a stable stochastic long-term growth rate (?_s = 1). 

Logging and its impact were documented at the Core 6 site that is near the larger Core 3 population found in the same regional coastal forest. We used data on the species identification and DBH of the logged host species to simulate the effects of two types of logging regimes on T. undulatum. Specifically, T. undulatum that were growing on species targeted by loggers were marked as dead, no matter the DBH of the host plants. The targeted species are Bucida buceras, Tabebuia angustata, and Annona glabra. A second matrix was created that used both host species and a minimum DBH measure from the nearby field observations to attribute mortality to the T. undulatum. We used the following minimum DBH sizes per logged host species: B. buceras = 18 cm, T. angustata = 22 cm, A. glabra = 23 cm. The ?_s was calculated for both logging scenarios to determine the frequency of logging under both regimes that would allow projected long-term persistence (?_s = 1) (Appendix S1: Table S2).

Citations

• Borrero, H., J. C. Alvarez, R. O. Prieto, and H. Liu. 2022. Comparisons of habitat types and host tree species across a threatened Caribbean orchid’s core and edge distribution. Journal of Tropical Ecology 38: 134–150.
• Caswell, H. 2001. Matrix population models: Construction, analysis, and interpretation (2nd ed.). Sinauer Associates. 
• Cochran, M. E. and S. Ellner, S. 1992. Simple Methods for Calculating Age‐Based Life History   Parameters for Stage‐Structured Populations: Ecological Archives M062-002. Ecological monographs 62 :345–364.
• Landsea, C. W. and J. L. Franklin. 2013. Atlantic hurricane database uncertainty and presentation of a new database format. Monthly Weather Review 141: 3576–3592.
• Lefkovitch, L. P. 1965. The study of population growth in organisms grouped by stages. Biometrics 21: 1-18.
• Stott, I., S. Townley, and D. J. Hodgson. 2011. A framework for studying transient dynamics of population projection matrix models. Ecology Letters 14: 959–970.
• Stott, I., D. J. Hodgson, and S. Townley. 2012. popdemo: an R package for population demography using projection matrix analysis. Methods in Ecology and Evolution 3: 797–802.
• Stubben, C., and B. Milligan. 2007. Estimating and analyzing demographic models using the popbio package in R. Journal of Statistical Software 22: 1–23.

The following text files were imported and read by R 2.4.4 Software Package using the code script from the R file "Trichocentrum_undulatum_HB"

The text file titled "tricho8yrs_feb25" used for "Part I" and "Part IV" labelled sections of the R script. "Part I" created matrices using the demographic data for each site, census, and plant identification found in the text file. "Part IV" was the transient analysis that projected the short-term dynamics and used the most recent field observations or stage class vector for each population.

The text file "tricho_matrices_nov29" was used for the stochastic simulation of hurricane introduction, logging introduction, episodic recruitment introduction, and scale insect mortality removal at select censuses/populations. The file was used in the "Part II" and "Part III" labelled section in the R script. 

The following text files were imported and read by R 2.4.4 Software Package using the code script from R files "tricho.mark.hur" and "Tricho_markov_hurriane"

The text file "Tricho_markov_hurriane" contains the state of environments for the hurricane probability scenarios.

The text file "tricho.mark.hur" contains the matrices that will be pulled for each phase as described in the above methods "four major transition phases gathered from the Peripheral population censuses: (1) phase I, the hurricane year (census 5); (2) phase II, first year post-hurricane (census 6); (3) phase III, second and third year post hurricane (7 and 8); (4) phase IV, non-hurricane affected years (census 1, 2, 3, and 4) (Table 2; Appendix S1: Table S1)."