DOI: 10.5061/dryad.vhhmgqnxd
title: Demographic study of a tropical epiphytic orchid with stochastic simulations of hurricanes, herbivory, episodic recruitment, and logging

description: 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.
description: <strong>Field methods</strong>
 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 <em>T</em>.
 <em>undulatum</em> presented in this study, individuals were
 studied within ~1-km strips where <em>T</em>.
 <em>undulatum</em> occurrence was deemed representative of the
 site, mostly occurring along informal forest trails. Once an individual of
 <em>T</em>. <em>undulatum</em> 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
 <em>T</em>. <em>undulatum</em>, 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 <em>M</em>.
 <em>miamensis</em>, and whether an inflorescence was partially
 or completely affected by the fly, and whether there was other herbivory,
 such as <em>D</em>. <em>boisduvalii</em> 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 <em>M</em>. <em>miamensis</em>. 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. <strong>Matrix
 modeling</strong> <em>Definition of
 size-stage classes</em> 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.
 <em>Matrix elements and population vital rates
 calculations</em> 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 &lt; 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. <em>Deterministic
 modeling analysis</em> 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 (?<sub>s</sub>), 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). <em>Stochastic
 simulation</em> 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. <em>Trichocentrum</em>
 <em>undulatum</em> 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 ?<sub>s</sub> = 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 <em>D</em>.
 <em>boisduvalii</em> 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 <em>D</em>.
 <em>boisduvalii</em> 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
 ?<sub>s</sub>. 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 (?<sub>s</sub>
 = 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 <em>T</em>. <em>undulatum</em>.
 Specifically, <em>T</em>. <em>undulatum</em> that
 were growing on species targeted by loggers were marked as dead, no matter
 the DBH of the host plants. The targeted species are
 <em>Bucida</em> <em>buceras</em>,
 <em>Tabebuia</em> <em>angustata</em>, and
 <em>Annona</em> <em>glabra</em>. 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
 <em>T</em>. <em>undulatum</em>. We used the
 following minimum DBH sizes per logged host species:
 <em>B</em>. <em>buceras</em> = 18 cm,
 <em>T</em>. <em>angustata =</em> 22 cm,
 <em>A</em>. <em>glabra</em> = 23 cm. The
 ?<sub>s</sub> was calculated for both logging scenarios to
 determine the frequency of logging under both regimes that would allow
 projected long-term persistence (?<sub>s</sub> = 1) (Appendix
 S1: Table S2).
 <strong>Citations</strong> 
 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
 <strong>38</strong>: 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 <strong>62 </strong>:345–364.
 Landsea, C. W. and J. L. Franklin. 2013. Atlantic hurricane
 database uncertainty and presentation of a new database format. Monthly
 Weather Review <strong>141</strong>: 3576–3592.
 Lefkovitch, L. P. 1965. The study of population growth in
 organisms grouped by stages. Biometrics <strong>21</strong>:
 1-18. Stott, I., S. Townley, and D. J. Hodgson. 2011.
 A framework for studying transient dynamics of population projection
 matrix models. Ecology Letters <strong>14</strong>:
 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
 <strong>3</strong>: 797–802. Stubben, C.,
 and B. Milligan. 2007. Estimating and analyzing demographic models using
 the popbio package in R. Journal of Statistical Software
 <strong>22</strong>: 1–23.
description: 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)."