Bayesian Estimates of Transition Probabilities in Seven Small Lithophytic Orchid Populations: Maximizing Data Availability from Many Small Samples

Predicting population dynamics for rare species is of paramount importance in order to evaluate the likelihood of extinction and planning conservation strategies. However, evaluating and predicting population viability can be hindered from a lack of data. Rare species frequently have small populations, so estimates of vital rates are often very uncertain due to lack of data. We evaluated the vital rates of seven small populations from two watersheds with varying light environment of a common epiphytic orchid using Bayesian methods of parameter estimation. From the Lefkovitch matrices we predicted the deterministic population growth rates, elasticities, stable stage distributions and the credible intervals of the statistics. Populations were surveyed on a monthly basis between 18–34 months. In some of the populations few or no transitions in some of the vital rates were observed throughout the sampling period, however, we were able to predict the most likely vital rates using a Bayesian model that incorporated the transitions rates from the other populations. Asymptotic population growth rate varied among the seven orchid populations. There was little difference in population growth rate among watersheds even though it was expected because of physical differences as a result of differing canopy cover and watershed width. Elasticity analyses of Lepanthes rupestris suggest that growth rate is more sensitive to survival followed by growth, shrinking and the reproductive rates. The Bayesian approach helped to estimate transition probabilities that were uncommon or variable in some populations. Moreover, it increased the precision of the parameter estimates as compared to traditional approaches.


Introduction
A large amount of effort has been placed in understanding population dynamics of orchids [1][2][3][4][5][6][7][8], focusing mainly on conservation and management decisions [9][10][11][12][13]. Almost all studies are of terrestrial temperate species, whereas approximately 70% of orchid species are epiphytic. Insufficient information on the basic biology and dynamics of the species often limits attempts to conserve small populations. When population sizes are small, decisions frequently rely on intuition or data from other populations or species. Parameter estimates are usually uncertain when data are limited. In some cases, it is not possible to even construct meaningful credible intervals [14]. Traditional approaches using standard statistical methods for evaluating population parameters can be limiting [14][15][16][17]. Here we use Bayesian methods to estimate the parameters of an orchid population model by using information from a number of populations simultaneously.
One value of the Bayesian approach is that it takes advantage of a priori information and considers the probability distribution for the parameter of interest [18]. Moreover, the 95% Bayesian credible interval for a parameter can be interpreted as intended by most biologists as the interval with a 0.95 probability of containing the true value [19][20]. A further advantage of the Bayesian approach is that the parameters of arbitrarily complex models can be estimated with relative ease [17][18], [21].
As part of a long-term objective of evaluating alternative approaches for the management of small hyper-dispersed species in lithophytic and epiphytic habitats, we evaluated the transition probabilities of a matrix population model seven riparian populations of a lithophytic tropical orchid from the Caribbean. Specifically we asked 1) what are the transition probabilities among life cycle stages, 2) which of the transitions have the highest and lowest elasticities, 3) what is the expected stable stage distribution, 4) how different are populations from a stable stage distribution, and 5) how different are populations from different riparian environments? In answering these questions, we focused on the size and precision of the parameter estimates.

Model system
The model organism Lepanthes rupestris Stimson is an epiphytic and mostly lithophytic (on boulders) endemic to Puerto Rico, and found in the Yunque National Forest (YNF). No specific permits were required for research as this species is neither endangered or protected (no removal or manipulation of plants were performed). For further information contact, El Yunque National Forest, HC-01, Box 13490, Rio Grande, PR, 00745-9625, USA. The study was carried out along two rivers, Quebrada Sonadora (latitude 18u 199170 N, longitude 65u 499 110 W) and Quebrada Grande (latitude 18u 199 090 N, longitude 65u 499 300W). These rivers are first order tributaries of the Espíritu Santo River in north-eastern Puerto Rico. The sections studied range between 400 and 500 m in elevation and occur in secondary mature 'Tabonuco Forest', dominated by Dacryodes excelsa [22].
For the purposes of this study, a population is defined as the set of individuals growing on a single boulder (lithophyte) or tree (epiphyte). For ease of communication we will use the term phorophyte to refer to both (tree and boulder) host types. There is a mean of 45.368.1 (s.e.) adults per phorophyte and a mean distance of 4.8 m61.3 (s.e.) between phorophytes [23]. Most populations have fewer individuals than the mean (Median 26); [6], [23]. Populations are clumped in space, small and separated from each other by short distances. In general, populations on separate phorophytes are genetically distinct and gene flow among populations appears to be limited [24] suggesting that each population has limited interactions with its neighbors. Lepanthes rupestris is self-incompatible and is protandrous [25]. Although much effort has been invested to discover the mechanism for transferring the pollinaria, pollinators are unknown for this species. The only known pollinator for one of .800 species of the Neotropical genus is described as a black winged fungus gnat [26] and suggests that pollination occurs through pseudocopulation by male black winged fungus gnats.

Description of Field Site
Average annual precipitation is 3460 mm but rainfall is seasonal, with higher values from May through December relative to the January to April period [27]. Water discharge in these streams is highly variable and closely follows precipitation events [28]. Large storms increase water discharge to levels that are much higher than normal and are capable of moving boulders and carrying large pieces of debris downstream [28]. The sections of Quebrada Sonadora are much wider (approx. 10 m) than those of Quebrada Grande (approx. 3 m). Most orchid populations are in the understory in the latter river and receive only indirect sunlight compared to populations along Quebrada Sonadora [29]. The size of plants is directly correlated with light availability in ex situ and in situ experiments and influences directly the production of flowers and consequently male and female reproductive success [29][30]. However, plant growth was non-linearly related to light availability with maximum growth occurring in the range 1 to 5 mol m 22 day 21 [29].

Data collection
A total of 1346 individuals of L. rupestris in seven populations were tagged and monitored at monthly intervals for 18 months or more throughout the sampling period in two river systems, Quebrada Sonadora and Quebrada Grande. The number of individuals varied among populations (Table 1). Analyses were performed using demographically structured models in which between life-history stages between successive census dates (monthly surveys). G = probability of passing from one stage to another, P = probability of surviving the time period and remaining in the same stage. Arrows labelled with G ij 's indicate movement of an individual from one stage to another. P ij 's indicate individuals remaining in the same stage at the next census date. The arrow labelled F ij indicates the production of new individuals by sexual reproduction by a specific stage. S = seedlings, J = juveniles, A 0 is a non-reproducing adult stage, representing individuals that have reproduced at some point in their life but which are presently not reproducing, A 1 = small reproductive adults (only one active inflorescence), A 2 = large reproductive adults (two or more active inflorescence). New individuals are introduced in the populations only from reproducing individuals. doi:10.1371/journal.pone.0102859.g001 individuals were grouped into distinct classes. Readers unfamiliar with the approach can read Caswell [31] or Akçakaya et al. [32] for a brief explanation of the methods and assumptions. Lepanthes rupestris is a small caespitose species of approximately to 15 cm tall; individuals have mainly 2 to 6 slender stems with a solitary coriaceous leaf [33]. Inflorescences are produced and are active throughout the year. Flower production on an inflorescence is sequential from the base, each flower last approximately one and half week (dependent of the environment) and fruits (globose ca. 4 mm; with thousands of seeds) are persistent for 6 weeks. At each survey, individuals were recorded as being in one of five stages, defined as follows ( Figure 1): (i) seedlings: small plants without petioles on any leaf, (ii) juveniles: individuals with at least one lepanthiform sheath on the petiole and no current or previouslyproduced inflorescences (these are persistent) (iii) non-reproducing adults: individuals that were not currently flowering, but carried dried inflorescences from a previous flowering event, (iv) small reproductive adults: individuals with one photosynthesizing (green) inflorescence with or without open flowers or fruits, (v) large reproductive adults: individuals with two or more reproductive shoots. Reproductive effort through female success (fruit set) by all reproducing individuals was noted at every recording date. Fruits remain on plants for about 1.5 months, as each fruit was monitored individually by identifying the leaf, inflorescence and position fruit on the inflorescence production was not overestimated. Consequently, it was possible to account for all fruits produced during the sampled period.

Statistical Analysis
Data on monthly survival and transitions among the adult age classes were analyzed with multinomial models. For example, nonreproductive adults could remain as non-reproductive adults (n), or transition to small (s) or large reproductive adults (l), or die (d). Mortality of adults in stage i, month j and population k (m k,j,i ) was modelled using logistic regression with categorical variables for stage, month and population. The regression took the form logit(m k,j,i )~stage i zmonth q½k,j zpop k , where stage i is an effect describing how mortality varies among the three adult stages, month q[k],j is an effect describing monthly variation in mortality in population k that depends on which stream it is in (q[k] took values of 1 or 2), and pop k is an effect  describing how mortality varies among the different populations.
The variable month q[k],j was treated as a random effect while the other three variables were treated as fixed effects. Within the two streams, monthly variation in mortality was assumed to be the same for all populations (i.e., within a given month j, populations in the same river shared the identical value for the random effect month q[k],j ), while no such correlation was assumed among streams.
Conditional on survival, we assumed that adults reproduced with a particular probability, and that reproductive adults were large or small with a particular probability. These probabilities varied among the type of adult (e.g., adults that were nonreproductive in one month and those that were large reproductive adults had different probabilities of reproducing in the next month) and among populations. However, we assumed that there was not monthly variation in these transitions (conditional on survival). Thus, transitions for an adult of class i (taking values n, s, l or d, that specify non-reproductive, small reproductive, large reproductive or dead individuals, respectively) to class n, s, l or d in month j and population k were specified by the equations where r k,i is the probability that adults of class i in population k are reproductive in the next month given that they survive, and s k,i is the probability that adults of class i in population k are large in the next month given that they survive and become reproductive. We specified, not surprisingly, that dead individuals (i = d) remained dead (i.e., p k,j,i,d = 0 for i = n, s and l, and p k,j,d,d = 1), but estimated the other transition probabilities using the available data.
Data on transitions between adults and the younger stage classes (juveniles and seedlings) and among the younger classes were available on approximately 13-monthly time scales. Conversely, juveniles had sufficient time to develop into small or large reproductive adults, but may remain as juveniles or transition to non-reproductive adults. These assumptions provided realistic  limits on the rates of development of individual plants from seedlings to adulthood. Reproduction rates were estimated using data on occurrence of new seedlings and juveniles in the populations. The data indicated that large reproductive adults produced 7.8 times as many fruits as small reproductive adults and consequently we assumed the ratio was in concordance with the production of seedlings and juveniles. We assumed that this ratio was the same in all populations. Similarly, the ratio of seedlings to juveniles being produced was assumed to be the same in all populations. The actual rate of seedling and juvenile production by reproductive adults was permitted to vary among the seven populations. We assumed that seedlings and juveniles that became non-reproductive adults in a 13-month period (i.e., had evidence of having produced fruits) had contributed to production of seedlings and juveniles at the same rate as small adults. Over this time period, we assumed that seedlings could become small reproductive adults, but did not have sufficient time to become large reproductive adults.
The survival of seedlings and juveniles was modeled using logistic regression, with the rate varying among populations, and between these two stages. Thus, survival of seedlings in population k was given by logit(s y,k ) = a k , and survival of juveniles was given by logit(s j,k ) = a k +b.
Conditional on the seedlings and juveniles surviving, the probabilities of progressing to older stages were expressed as a series of conditional probabilities. Thus, the probability of a seedling (class y) becoming a juvenile (class j) in population k was equal to p k,y,j = s y,k g y t s,j , where g y is the probability of the seedling progressing to an older stage conditional on surviving, and t s,j is the probability of a seedling becoming a juvenile conditional on the seedling surviving and not remaining a seedling.
For seedlings, the probability of becoming a non-reproductive adult was p k,y,n = s y,k g y (1-t s,j )t s,n , where t s,n is the probability of a seedling becoming a non-reproductive adult conditional on it surviving. Finally, the probability of a seedling becoming a small reproductive adult was p k,y,n = s y,k g y (1-t s,j )(1-t s,n ). Similar expressions were used for the transition probabilities of juveniles to the juvenile, non-reproductive, small reproductive and large reproductive stages.
The parameters of the model were estimated using the observed data on transitions for the seven populations. The estimation procedure was conducted in WinBUGS, a program for conducting Bayesian analyses with Markov chain Monte Carlo methods [37]. This was used because it permitted a relatively easy analysis of the complex hierarchical model [17] and accounted for instances when data were missing. Missing data such as the absence of detection of transition from one stage to another can be estimated for a specific population from the whole data. In addition when events are rare in a population or sample size are small, parameter estimates can be unrealistic if estimated from the data for the single population, hierarchical model uses the whole data set to calculate the best parameter estimates for each specific populations.  Consequently, Winbugs allows for complex hierarchical models, which in this case estimates transition probabilities from logistic regression for survival and transitions. Thus one advantage is that it estimates simultaneously the parameters of survival and transitions for all populations while traditional methods estimate these individually for each population [31].
To reflect a lack of prior information, vague prior distributions were used for the parameters. Uniform distributions between zero and one were used for probabilities (e.g., g y , t y,j , etc), and normal distributions with means of zero and standard deviations of 1000 were used for the regression coefficients (e.g., a k , b, stage i and pop k ,). The large standard deviation for the prior meant that it had essentially no influence on the parameter estimate. The random effects for monthly variation in adult mortality (month j ) were modeled as deviates drawn from a normal distribution with a mean of zero and standard deviation that had a wide uniform prior distribution.
Code for conducting the analysis is given in the Appendix. A total of 10000 samples from the posterior distributions of the parameters were taken after discarding the initial 5000 samples as a burn in. The sampled parameter values were exported from WinBUGS and analyzed in the program Mathematica to calculate the asymptotic growth rate (dominant eigenvalue of the 13monthly transition matrix), stable stage distribution and elasticities of the parameter values as defined by Caswell [31]. These population parameters were calculated within Mathematica for each set of parameter values generated from WinBUGS, and then their mean and standard deviation were calculated to characterize their posterior distributions. Monthly survival rates of adults were converted to 13-monthly rates by raising a sub-matrix of monthly transition probabilities to the power of 13 (i.e., we did not allow multiple reproduction events by adults within a single 13-month period). Our model does not contain density dependence, catastrophes and genetic effects.
For those unfamiliar with the Bayesian approach we recommend that they explore the terminology and methodology with these approachable books [21], [38][39][40].

Variation in reproductive potential
Flower production varied greatly among individuals and populations ( Table 2). In a few populations most adult individuals produced flowers (pop 5 and 7) while in other populations the number of individuals that failed to produce flowers was significantly higher (pop 1). Pollinaria removal (evidence of pollinator activity) and fruit set showed a similar pattern of large variation among individuals and populations. This suggests that pollinator activity is different among populations.
In all populations larger plants had a higher probability of bearing more fruits than small reproductive adults. A small reproductive plant had a 1 to 4% probability of bearing a fruit in any of the months sampled, while a plant having two or more active inflorescences had a 7 to 20% probability of bearing fruits.  In spite of this, most adult plants when reproductive were small, only a very small fraction of reproductive adults had two or more inflorescences (see below, section on stable stage distribution).

Transition Probabilities
The stage with the highest probability of death during the survey period was dominated by seedlings and was consistent for all populations (10-51%; Table 3 to Table 9). The juveniles and non-reproductive adult stages shared the second highest probability of death, while large reproductive adults were least likely to die (3-20%). The transition probabilities varied among the seven populations. For example, the likelihood of a seedling staying in the same stage in the next time period ranged from 19 to 34%.
The stage that appears to impede growth in the life history of Lepanthes rupestris is the juvenile stage where, in most of the populations, individuals remained juveniles in the next survey period (48-72%). Large adults did not stay in that stage for long; only in populations one and five was the likelihood of remaining as a large adult greater than 50% while in population three the likelihood of remaining in that stage was much lower (31%). Nonreproductive adults had a high probability of becoming reproductive (flower production) in the next time period in almost all populations except population three, where non-reproductive adults often remained in that stage (Table 3 to Table 9).

Asymptotic population growth rates of individual populations
The asymptotic population growth rates of Lepanthes rupestris had a range of 0.98 to 1.01 (Figure 2). In three of the populations the asymptotic population growth rate cannot be distinguished from stability (1, 4 and 5) because the 95% credible interval overlapped 1.00, while two populations had asymptotic growth rate below one by approximately 1-2% (5 and 6) and two populations had asymptotic growth rates above one (2 and 3). The positive growth in population 2 is small and suggests a maximum increase of less then 1%, while population 3 was increasing by about 1-2%.
Sites on both river systems (Quebrada Grande, pop 1-3; Quebrada Sonadora, pop 4-7) had similar population growth rates (Fig. 2), although the highest population growth rate occurred on Quebrada Grande, the shadier river, and the two lowest population growth rates occurred on Quebrada Sonadora.

Elasticities
In all populations, modifying the parameters for recruitment (F ij ), survivorship (P ij ) and growth (G ij ) of seedlings would have little impact on population growth rates as most elasticities are small (elasticities #0.005; Table 10 to Table 16). In almost all cases the largest elasticities were observed for the proportion of large reproductive adults remaining reproductive (with a maximum elasticity of 0.276). In six of the populations, the largest elasticities are associated with the larger class sizes of adults, and mainly influenced by small and large reproductive adults remaining reproductive and growing to the larger size class. The exception is population three where the proportion of juveniles remaining juveniles has the largest elasticity, with low values for all other transitions. The sum of the elasticities of the small and large adults was comprised between 0.59 to 0.71, except for population three where it comprised a smaller fraction (0.46).

The expected stable stage distribution
If populations of L. rupestris do attain stability then the expected stable stage distribution of the different populations suggests that the juveniles (in four populations) or the large adults (in three populations) are likely to be the most common of the five stages (Table 17). In all populations, the frequency of seedlings is expected to be small (,1-4%). The expected proportion of reproductive adults (small and large adults) varies extensively among populations from 27-61%. In all populations of Quebrada Grande the juvenile stage comprised the largest proportion of individuals at the stable stage distribution while at Quebrada Sonadora the dominant stage varied among population.
The observed frequency distribution at the beginning of the survey was significantly different from the predicted stable stage distribution for all populations (all chi square p,0.0001), with the number of observed seedlings higher than expected. The rarity of populations observed at stable stage distribution in epiphytic systems and orchids in general suggest that there is little evidence that these orchids ever attain and sustain stable stage distribution [1][2][3][4][5][6][7], [24], [34], [36], [41]. It is possible that for many orchid species populations may always be in a transient phase [42][43][44].

Transition rates
In general a large proportion of seedlings stayed seedlings (19-34%) for more than a year and most juveniles did not growth to the adult stage and remained in the same stage (48-72%). Because of the differences among the width of the two rivers, Quebrada Sonadora and Quebrada Grande, a priori we expected life history differences between the two sites. We found that male and female fitness did vary among populations but with no obvious environmental differences [45] except that populations at Quebrada Grande were generally shadier, and likely more humid, than at Quebrada Sonadora, a wider river with more direct sunlight. If pollinators are fungus gnats as we suspect, then local environmental conditions such as humidity, quantity and quality of forest litter, and soil organic matter may control their  populations and abundance. Consequently we should see more pollinator activity at Quebrada Grande because of the shadier environment. There was little or no difference in pollinaria removal between the two river systems, but plants at Quebrada Grande experienced higher fruit set. Why did we observe this difference? The answer may rest more with flower response to humidity than any effects on pollinators. Flowers of Lepanthes are sensitive to dry conditions; low humidity reduces flower longevity (Tremblay unpublished data). In addition, stigmatic maturity depends on flower age (L. rupestris is protandrous), as does the probability of successful pollination and fruit set [25]. Consequently, flowers of Quebrada Sonadora populations may be shorter-lived with reduced opportunity for pollinations. This may explain the tendency for higher asymptotic population growth rates at Quebrada Grande. A constant observed in all populations surveyed is that the transition from juvenile to adulthood is appears to be one of the bottlenecks in the life history of the orchid as this transition was rarely observed during the survey period.
Out of a total of 4343 juveniles surveyed (sum of all populations and all months) only 88 of these became adults. What is the cause of this low transition rate to the adult stage but it is likely to be a consequence of multiple environmental (abiotic and biotic interactions) and stochastic factors. The transition between the three adult stages was frequent, and large reproductive individuals were not likely to be in that stage during the next survey period (36-69%). Reproductive effort in all of the L. rupestris populations was low and most often limited to a small fraction of adults. Reproductive success is pollinator limited and is highly skewed towards a few individuals [24]. Once the adult stage is attained, mortality is low.
In an extensive review of the population ecology of orchids mediated by environmental factors, Light et al. [7] list many studies. However, most of these are short-term studies and descriptive, which fails to use the information to predict population persistence. Most orchid population dynamic studies have been restricted to terrestrial species (for a review see) [46]. Comparative data on epiphytic and lithophytic species are limited to four species [34], [47][48][49][50][51][52].
Large individuals have high reproductive potential and lifespan, but these individuals are rare in the population. The lifetime reproductive success of Lepanthes is extremely skewed where most individuals have low (or zero) reproductive success and a few have very high lifetime reproductive success [24]. In an 11-year study of the related Lepanthes caritensis Tremblay and Ackerman, some adults may survive for more than 20 years but most individuals live less than 4 years [36]. In many species of orchids, larger individuals have higher reproductive potential or success [53]. There is a correlation between life stages and size of individuals in Lepanthes; generally leaf size and leaf number increase as an individual progresses from juveniles to small reproductive adults  and large adults. Thus a complicating life history variable is that a large proportion of reproductive orchids become smaller in the next survey period suggesting that flower production or fruit production may be costly. A cost of reproduction has been observed in many species of orchids [53]. In epiphytic species, subsequent plant size (pseudobulb or leaf number or size) is often smaller after a reproduction. In terrestrial species, plants are often dormant the next growth period or emerge as non-reproductive individuals [13], [54], but this is not always the case [55]. In general, survival and growth were the most important demographic processes in all of the populations. Stasis (sum of P ij ) had the largest of the elasticities in all seven populations (0.10-0.28). While relative contribution of fecundity to change in population growth was negligible for all populations (0.004-0.011). Elasticites of all populations are dominated by growth and survival with fecundity having little impact. This pattern of dominance of elasticities from growth and survival suggest that Lepanthes rupestris are more similar to iteroparous forest herbs and not from open habitats [56].
Comparative data from orchids are rare [13], [34], [52], [57]. However, most studies suggest that survivorship is important. For example, in the terrestrial Ophrys sphegodes over 87% of the elasticities is explained by survivorship of reproductive individuals [57]. In the epiphytic Lepanthes caritensis the highest elasticities are found in stasis of the non-reproductive adults (.0.52%; Tremblay 1997), while in the rare species Lepanthes eltoroensis the highest elasticities varied among the five populations. In three populations the survivorship and stasis of non-reproductive adults was the highest and in two populations the stasis of reproductive adults [34]. It appears that promoting flowering production and remaining a reproductive adult would have the largest impact on population growth, however, the relative importance of these cannot be determined a priori. In the terrestrial Prasophyllum correctum the most important demographic process is the dormant stage (.0.31) [13].

Bayesian advantage
Our paper illustrates one of the benefits of using of a Bayesian approach for estimating transition probabilities; relatively complex statistical models can be analyzed [17]. These analyses can use information from a wide range of sources simultaneously to fit complex hierarchical models. Because information from several populations is used concurrently, the resulting parameter estimates tend to be more realistic compared to the case where each population is analyzed individually using, for example, traditional methods as in Ebert [58] or Caswell [31]. An example of this is that estimating individual transitions for juveniles for two of the populations was hampered because we observed no deaths during our survey period (data not shown). This result is not unusual when sample sizes are small. Similarly in one of the populations we observed only 1 of 37 (3.7%) individuals that grew from a juvenile to an adult stage, while in another population 38 of 778 (1.3%), which are the best maximum likelihood of the observed data. The Bayesian estimation of this transition can be further refined by  taking advantage of all the data from all populations from gave estimates that range from 5.93% to 8.96%, and more importantly, credible intervals of the transition probabilities account for estimation errors. The differences in transition rates could easily be assigned to difference in sample size and not solely as a result of environmental or genetic differences among populations.