Fallen Pillars: The Past, Present, and Future Population Dynamics of a Rare, Specialist Coral-Algal Symbiosis

Rare and ecologically specialized species are at greater risk of extinction. Rarity in terms of low genotypic diversity may be obscured in sessile foundation species that can reproduce asexually. Asexual propagules are often only distinguishable from sexual recruits through molecular tools. Thus, molecular markers are necessary to assess genotypic variation and population structure in clonal organisms such as corals. The global decline of corals has necessitated marker development for improved conservation of rare coral species. We infer past demographic changes, describe modern population structure, and quantify asexual reproduction of the uncommon Caribbean pillar coral, Dendrogyra cylindrus and its endosymbiotic dinoflagellate, Symbiodinium dendrogyrum using de novo microsatellite markers. Results show that D. cylindrus comprises three distinct populations in the Caribbean whereas the symbiont was differentiated into four populations. Thus, barriers to gene flow differ between host and symbiont. In Florida, host and symbiont reproduced mainly asexually, yielding lower genotypic diversity than predicted from census size. Models of past demographic events revealed no evidence of historical changes in population size, consistent with the geological record of D. cylindrus indicating it has been rare for hundreds of thousands of years. The most recent global thermal stress event triggered a severe disease outbreak among D. cylindrus in Florida, resulting in a precipitous population decline. Projections indicate a high likelihood that this species may become locally extinct within the coming decades. The ecosystem consequences of losing rare coral species and their symbionts with increasingly frequent extreme warming events are not known but require urgent study.


21
Rare and ecologically specialized species are at greater risk of extinction. Rarity in terms of low 22 genotypic diversity may be obscured in sessile foundation species that can reproduce asexually. 23 Asexual propagules are often only distinguishable from sexual recruits through molecular tools. 24 Thus, molecular markers are necessary to assess genotypic variation and population structure in 25 clonal organisms such as corals. The global decline of corals has necessitated marker 26 development for improved conservation of rare coral species. We infer past demographic 27 changes, describe modern population structure, and quantify asexual reproduction of the 28 uncommon Caribbean pillar coral, Dendrogyra cylindrus and its endosymbiotic dinoflagellate, 29 Symbiodinium dendrogyrum using de novo microsatellite markers. Results show that D. Ecologically rare species are predicted to be more vulnerable to environmental change 44 and are at greater risks of extinction with shifts in climate (Caughley, 1994;Davies et al., 2004;45 McKinney, 1997). Species are rare because they may inhabit narrow geographic ranges, occupy 46 few specific habitats, and/or exhibit low abundance in nature (Rabinowitz, 1981). It is sometimes 47 presumed that rare species tend to be competitively inferior, although work on sparse and 48 common grasses does not support this conclusion (Rabinowitz et al., 1984). Regardless, there are 49 obvious consequences to having low population densities, including difficulties finding a mate 50 (Stephens & Sutherland, 1999) and vulnerability to genetic drift (Ellstrand & Elam, 1993). 51 However, there are many rare species that are adapted to low population densities (de Lange & 52 Norton, 2004;Flather & Sieg, 2007), and these species persist for long periods of evolutionary 53 time. 54 The persistence of rare species is further challenged when they are obligate partners of 55 specific symbionts. The endangered terrestrial orchid Caladenia huegelii associates with a 56 specific mycorrhizal fungus throughout its range (Swarts et al., 2010). This specificity between 57 partners has caused C. huegelii to be rare due to the limited suitable environmental conditions of 58 the mycorrhiza. Thus, the strict niche characteristics of one partner in a symbiosis may drive the 59 scarcity of the other. Further, intra-specific diversity in both partners and fidelity of genotype-60 genotype associations can play a role in how the symbiosis responds to changing conditions 61 (Baums et al., 2014b;Parkinson et al., 2015;Parkinson & Baums, 2014). 62 Because sessile organisms often reproduce by asexual fragmentation, the relative rarity of 63 these species may be obscured by an apparently large census size. Thus, molecular markers are 64 Caribbean, some species were consistently more abundant (e.g. Orbicella faveolata), while 88 others were rare, such as the pillar coral Dendrogyra cylindrus (Edmunds et al., 1990; S. Miller 89 et al., 2013;Steiner & Kerr, 2008;Ward et al., 2006). The fossil record supports the relative 90 rarity of D. cylindrus through time, although there is localized evidence that pillar corals were 91 more prevalent on Pleistocene reefs (Hunter & Jones, 1996 (Finney et al., 2010). Increasing rarity of one species will likely lead to 99 subsequent rarity in the other mutualistic partner. While it is unclear which partner is driving the 100 rarity of the D. cylindrus-S. dendrogyrum mutualism, assessing levels of within-species diversity 101 is imperative to understanding the functioning of and threats to this association. If co-dispersal 102 between partners is occurring, then the population structure of both symbiotic species will be 103 Table 1). However, after conducting preliminary analyses using this set of markers, it became 156 clear that additional markers would be necessary to increase our power of identifying clonal 157 strains. See the supplement for details on de novo primer design and amplification conditions. 158 The combined set of 15 microsatellite markers (Table 2) was used to assess clonal and 159 population structure of the algal symbiont. 160

Analysis of Population and Clonal Structure 161
Matching multilocus genotypes were identified using the Data Subset option in 162 GENALEX vers 6.5 (Peakall & Smouse, 2006 (Pritchard et al., 2000). After initial testing, the admixture model 176 with no location prior was used with correlated allele frequencies. The prior number of 177 populations (K) was set from 1 to 6 with five replicate runs per K, a burnin of 100,000 and 178 1,000,000 Markov Chain Monte Carlo repetitions after the burnin using the package 179 'PARALLELSTRUCTURE' (Besnier & Glover, 2013) in R v3.4 (R Development Core Team 2017). 180 The optimal value for K was identified using the delta K method (Evanno et al., 2005) 181 implemented in the online program STRUCTURE HARVESTER (Earl, 2012). CLUMPAK was used to 182 identify the consensus of inferred clusters for the different replicates of each K value and to 183 visualize the results (Kopelman et al., 2015). Additional estimators of the number of optimal 184 populations (posterior probability, MedMeaK, MaxMeaK, MedMedK, and MaxMedK) were 185 applied (Pritchard et al., 2009;Puechmaille, 2016). Because uneven sampling has been shown to 186 underestimate the true K (Puechmaille, 2016), larger populations were randomly subsampled to 187 create more even sample sizes across the five regions. Global F ST , or the proportion of genetic 188 variation between subpopulations (Wright, 1951), was calculated using GENODIVE (Meirmans &  Analysis of molecular variance (AMOVA) was employed to test hypotheses about 194 population structure of host and symbiont (Excoffier et al., 1992). Samples were grouped by 195 region (identified using STRUCTURE and the optimal K value) and location (Belize, Florida, Turks 196 and Caicos Islands, USVI, and Curaçao). The AMOVA option in GENALEX vers 6.5 was used to 197 partition genetic variation from a distance matrix between regions and locations, using 9999 198 permutations and assuming an Infinite Allele Model. 199 GPS coordinates for all Florida colonies were used to create a geographic distance 200 matrix, which was compared to the genetic distance matrix for D. cylindrus. Mantel tests of 201 isolation by distance and spatial autocorrelation analyses were conducted in GENALEX vers 6.5. 202 Because of the small number of samples, a finding of no population structure could be 203 due to a lack of power. Thus, simulations were conducted in POWSIM v.4.1 to assess the power 204 of the microsatellite data set to detect low levels of population differentiation. POWSIM 205 estimates the lowest level of differentiation (F ST ) that can be detected in a simulated population 206 with a minimum of 90% accuracy (Ryman & Palm, 2006). Significance was determined using 207 Fisher's exact test. 208

Population Demographic Modeling 209
The modeling software MIGRAINE v.0.5.2 (http://kimura.univ-210 montp2.fr/~rousset/Migraine.htm) was used to test for past changes in population size in the 211 Florida D. cylindrus population. Samples from other locations were not grouped with Florida 212 because population structure within a dataset yields inaccurate demographic parameter inference 213 (Leblois et al., 2014). Additional samples were included for this analysis to increase the power of 214 detecting historical population size changes (n=95 total). The null model assuming a constant 215 population size (OnePop) was run to estimate pGSM, which was close to 0.3. This parameter 216 describes the geometric distribution of mutation step size under a generalized stepwise mutation 217 model (GSM), which is most appropriate for microsatellite markers (Leblois et al., 2014). When 218 we ran the model assuming one past change in population size (OnePopVarSize), we fixed 219 pGSM at 0.3 and inferred the current (2Nµ) and ancestral (2N anc µ) population sizes, as well as 220 the duration of the contraction or expansion event (Dg/2N) and the ratio of current to ancestral 221 population size (N ratio ). We assumed a mutation rate (µ) of 5 x 10 -4 per locus per generation, a 222 value typical for microsatellite loci (Estoup & Angers, 1998). In addition, we also ran the 223 OnePopFounderFlush model, which assumes two past changes in population size and infers a 224 founder population size (2N founder µ) and two additional population size ratios (N cur N founder ratio and 225 N founder N anc ratio ). Significant population size changes were detected using 95% confidence 226 interval estimates of the population size ratios. If the value 1 (indicating the current and ancestral 227 population sizes were identical) was outside of the 95% confidence interval, then the size change 228 was significant. MIGRAINE uses sequential importance sampling algorithms (De Iorio & 229 Griffiths, 2004a, 2004b to estimate parameters of past demographic changes from population 230 genetic data. We used 2,000 points, 2,000 trees, and four iterations per run for the models with 231 population size changes. The null model was run using 2,000 points, 100 trees, and three 232 iterations. Two independent runs were conducted per model by changing the estimation seed. 233 Additional model settings are included in the supplement (Supplemental Table 2). 234

Projected Population Declines 235
Through the course of this study, we witnessed a severe decline in Florida D. cylindrus 236 (Neely et al. in prep). In Broward County, 86% of the known colonies were lost in two years 237 (Kabay, 2016). Severe thermal stress events such as this one are expected to occur annually by

Hardy-Weinberg Equilibrium and Linkage Disequilibrium 262
The 11 de novo microsatellite markers developed for D. cylindrus yielded few deviations 263 from Hardy-Weinberg equilibrium (Supplemental Table 3). None of the markers significantly 264 deviated from Hardy-Weinberg in more than one location. Few locus pairs showed significant 265 linkage disequilibrium, and none of the significant deviations was found in all locations. 266 Deviations from Hardy-Weinberg equilibrium could not be tested for the symbiont markers, 267 since S. dendrogyrum is haploid. Tests for linkage disequilibrium between the S. dendrogyrum 268 markers yielded many apparently linked locus pairs. However, this is expected in a species that 269 primarily reproduces asexually within its host, and thus recombination may be rare. Most paired 270 combinations of the symbiont loci were in equilibrium in at least one location, indicating that 271 these markers are not physically linked. Other locus pairs showed significant linkage 272 disequilibrium in just one location, without enough information available in the other locations to 273 discount linkage equilibrium. 274

Genetic Diversity of de novo Dendrogyra cylindrus microsatellite markers 275
The microsatellite markers developed for Dendrogyra cylindrus ranged in allelic 276 diversity from 6 to 15 and in effective allelic diversity from 2.047 to 6.849 (Table 3). Observed 277 heterozygosity levels ranged from 0.399 to 0.954, expected heterozygosity within subpopulations 278 ranged from 0.533 to 0.888, and total heterozygosity ranged from 0.686 to 0.895 (Table 3). The 279 total heterozygosity adjusted for sampling a limited number of populations ranged from 0.724 to 280 0.897 (Table 3). The heterozygosity values for all markers were relatively high (closer to the 281 maximum value of 1), indicating high genetic variability for these markers (Table 3). Inbreeding 282 coefficients for the 11 markers were between -0.19 and 0.361, with most values close to zero 283 (Table 3). The new markers differed in their genetic variability, with D15 showing the lowest 284 differentiation across all four measures (Table 3). D716 and D48 showed the most variability. 285 Overall, the markers were polymorphic enough to assess clonal and population structure in D. for genotypic diversity were low in both host and symbiont (albeit, slightly higher in the 312 symbiont), demonstrating that both species reproduced mostly asexually within a site. Genotypic 313 evenness was higher in the symbiont relative to the host, meaning that the host population was 314 dominated by relatively few highly replicated genotypes. 315 For the coral host, all of the clonal ramets of the same genet were contained within a 316 single collection site. This indicates limited dispersal abilities of asexual fragments of D. 317 cylindrus. A single coral genet dominated most of our sites, however, three sites had two genets 318 and one site contained three different genets ( Figure 4A). 319 A test of isolation by distance for D. cylindrus along the Florida Reef Tract (n=50) 320 revealed no significant correlation between genetic and geographic distance (r 2 = 0.0001). Spatial 321 autocorrelation analysis using the complete dataset of D. cylindrus genotypes (including clones) 322 in Florida (n=180) revealed significant positive spatial autocorrelation up to distances of 60 323 meters (Supplemental Figure 2A). However, when clones were removed, there was no 324 significant spatial autocorrelation (Supplemental Figure 2B). These results indicate that the cause 325 of the correlation between genetic distance and geographic distance over small distance classes is Florida are indeed the same species (S. dendrogyrum). Because such a small sample of symbiont 397 strains was included from both Belize and the Turks and Caicos Islands, it is possible that the 398 cluster containing these samples and some Florida samples is an artifact of STRUCTURE 399 (Puechmaille, 2016). 400

Historical Changes in Pillar Coral Population Sizes 401
Demographic modeling results revealed no evidence for past changes in pillar coral 402 population size. In both the one continuous population size change (OnePopVarSize) and the two 403 population size changes (OnePopFounderFlush) models, the upper limit of the 95% confidence 404 intervals for the ancestral population size (2N anc µ, both models) and the founder population size 405 (2N founder µ, OnePopFounderFlush model only) were undefined. The undefined upper limit meant 406 that many possible population sizes for the ancestral and founder population were equally likely, 407 including population sizes that were both larger and smaller than the 95% confidence interval for 408 the current population size (2Nµ). This was further evidenced by the broad profile likelihood 409 ratios for the ancestral and founder population sizes (Supplemental Figure 3), and the population 410 size ratios that included the value 1, indicating no significant change. To decide which of the 411 three models best described our data, we used the Akaike Information Criterion (AIC) calculated 412 based on the log likelihood and the number of parameters in each model (Anderson, 2007 However, the differences in AIC between models were very small, and thus it is possible that all 417 three models do not accurately describe reality. Still, the results from the models assuming past 418 demographic events failed to detect significant changes in population size, thus supporting our 419 conclusion of no past changes in D. cylindrus populations. 420

Projected Population Declines 421
The three projected rates of decline yielded a range of times to extinction for D. Over short distances (<60 m), Florida Dendrogyra cylindrus colonies were highly clonal 447 and showed positive spatial autocorrelation, indicating that the primary mode of reproduction 448 over this scale is asexual fragmentation (Supplemental Fig 2, Figure 3). These clonal structure 449 results are similar to those of the elkhorn coral, Acropora palmata (Baums et al., 2006)  compatible colonies is too low for successful sexual reproduction and population growth 466 (Courchamp et al., 1999). This study was limited to assess clonal structure in Florida alone, and 467 thus future work should measure genotypic richness, diversity, and evenness in other locations. 468 The Florida population of S. dendrogyrum is also clonal, with the same symbiont strain 469 often (but not always) found in all of the ramets of a coral genet. High fidelity between 470 individual genotypes of symbiotic partners has been found previously in the Caribbean elkhorn 471 coral (Acropora palmata) and its algal symbiont (Symbiodinium 'fitti') (Baums et al., 2014b). 472 While coral genets were restricted to one site each, five pairs of sites did share symbiont strains. 473 The sites that shared symbiont strains were close together ( Figure 5). Thus, the symbiont may be 474 able to disperse asexually over larger distances than its coral host, although successful long-475 distance dispersal appears to be rare. A typical fragment dispersal distance for D. cylindrus is 476 less than 60 meters, as evidenced by the spatial autocorrelation analysis (see below and 477 Supplemental Figure 2). 478 Similarly to the situation in D. cylindrus and S. dendrogyrum, the asexual dispersal 479 ability of Symbiodinium fitti was found to be significantly higher than its coral host, A. palmata fossil record (see Introduction). Abundance can be high in localized areas, however, due to 533 asexual fragmentation. This suggests that historically, D. cylindrus did not experience a range-534 wide decline since its first appearance in the fossil record during the late Pliocene/early 535 Pleistocene (Budd, 2000). 536 Through the course of this study, we witnessed a severe decline in Florida D. cylindrus 537 reef-building coral species in the Caribbean that will face local if not regional extinction. Global 563 change will continue to alter the composition of reefs if green house gas emissions are not 564 curbed. 565 The Consequences of Rare Species Loss 566 What are the consequences of losing rare species? Experimental removal of rare species 567 was shown to reduce ecosystem resistance to invasion by an exotic grass (Lyons & Schwartz, 568 2001). It is possible that less common species significantly contribute to the proper maintenance 569 of ecosystem function. Rare plant species were shown to more significantly impact nutrient 570 cycling and retention in an alpine meadow compared to more abundant species (Theodose et al., 571 1996). In a comprehensive study of species occurrence datasets from coral reefs, alpine 572 meadows, and tropical forests, rare species were repeatedly shown to predominantly support 573 vulnerable functions. These vulnerable functions were defined as ecosystem roles with low 574 redundancy, in that uncommon species with distinct trait combinations bolstered these particular 575 functions (Mouillot et al., 2013). Despite the high diversity in these ecosystems, abundant 576 species did not insure against the services lost by removing rare species. Thus, maintaining these 577 uncommon species is essential to overall community functional diversity. 578 When two rare species are combined in an obligate symbiosis, then the loss of one 579 species would yield coextinction of the other (Koh et al., 2004). It is possible that D. cylindrus 580 sexual recruits can survive by associating with S. meandrinium in the absence of S. dendrogyrum 581 but the latter has never been found in another coral host species (C. Lewis (Nei, 1987). §Hedrick's G' st is standardized relative to the maximum level of differentiation based on 953 the heterozygosity within subpopulations (Hedrick, 2005). ¶Contrastingly, Jost's measure of population 954 differentiation is independent of the within subpopulation diversity (Jost, 2008).