Introduction
Ancient lakes, continental bodies of water that persist on the landscape through multiple glacial cycles, often have a diverse endemic biota and therefore are ideal systems for studying evolutionary processes and their relationship to environmental change (Brooks Reference Brooks1950; Martens Reference Martens1997; Cohen Reference Cohen2003). Well-documented endemic faunas, such as the cichlid fish of Lake Tanganyika (Koblmuller et al. Reference Koblmuller, Sefc and Sturmbauer2008) and the invertebrate amphipod radiations in Lake Baikal (MacDonald et al. Reference MacDonald, Yampolsky and Duffy2005), have been intensely studied to better understand what produces biodiversity in ancient lake systems. Long-term stability (longevity) or periodic instability (landscape or climate change) are two alternative explanations for rich endemic diversity. Here, we investigate the morphology of the dominant phytoplankton species in the modern Lake Titicaca, the endemic diatom Cyclostephanos andinus, and its closely related morphospecies (hereafter referred to as the C. andinus complex), to determine whether there is evidence that the environmental history of the lake basin influenced the evolutionary trajectory of the species complex. To better characterize and understand the changes in the morphology of the C. andinus species complex in Lake Titicaca, we determined whether or not there are multiple discernible morphotypes of C. andinus, and we tested whether change in the body size of C. andinus exhibited stasis, a random walk, directional change, or punctuated change.
Diatoms are well preserved in the fossil record because of their silica cell wall. This attribute, coupled with their utility in reconstructing past environmental conditions, has led to the identification of at least 30,000, but possibly over 100,000, diatom species, largely on the basis of morphological characters (Mann Reference Mann1999; Mann and Vanormelingen Reference Mann and Vanormelingen2013). More recently, molecular studies have been used to complement morphological descriptions and sometimes have revealed cryptic or semicryptic species—those species that are morphologically similar but genetically distinct (Beszteri et al. Reference Beszteri, John and Medlin2007; Poulícková et al. Reference Poulícková, Veselá, Neustupa and Skaloud2010; Lundholm et al. Reference Lundholm, Bates, Baugh, Bill, Connell, Leger and Trainer2012). Species complexes of closely related taxa that are morphologically similar or indistinguishable often comprise a substantive component of the biodiversity in ancient lakes (Edlund et al. Reference Edlund, Levkov, Soninkhishig, Krstic and Nakov2006; Schön et al. Reference Schön, Poux, Verheyen and Martens2014; Wagner et al. Reference Wagner, Wilke, Krastel, Zanchettta, Sulpizio, Reicherter, Leng, Grazhdani, Trajanovski, Francke, Lindhorst, Levkov, Cvetkoska, Reed, Zhang, Lacey, Wonik, Baumgarten and Vogel2014). For those groups that leave a fossil record, and in systems where sedimentation is continuous over long periods, analysis of compositional, morphological, and genetic variation can be used to test evolutionary concepts over time intervals that bridge microevolutionary and macroevolutionary scales.
An ongoing debate in evolutionary biology is over the relative importance of phyletic gradualism versus punctuated equilibrium (Benton and Pearson Reference Benton and Pearson2001). Punctuated equilibrium (Eldredge and Gould Reference Eldredge and Gould1972; Gould and Eldredge Reference Gould and Eldredge1977) is a pattern defined by intervals of long-term stasis accompanied by short bursts of rapid evolution, in contrast to the incremental change characteristic of phyletic gradualism. Commonly, qualitative visual inspection of data has been used to determine whether phyletic gradualism or punctuated equilibrium is the dominant mode of change in an evolutionary lineage, and in some cases, opposing conclusions have been derived from the same data set, as exemplified by the interpretation of evolution in the fossils of Hyopsodus (Gingerich Reference Gingerich1976; Gould and Eldredge Reference Gould and Eldredge1977; Hunt Reference Hunt2008a). To alleviate the ambiguities in interpretation of evolutionary mode, Hunt (Reference Hunt2006, Reference Hunt2008b) proposed a likelihood framework to test competing models of phyletic evolution of fossil lineages, including directional evolution, random walks, stasis, and punctuated change.
The modeling of directionality in phyletic evolution, or what Hunt (Reference Hunt2006, Reference Hunt2008b) refers to as a “general random walk,” is a special case of a random walk, in which the mean has a nonzero step size (direction) for a trait over time. Conversely, an “unbiased random walk” is when the mean step size of a trait is equal to zero, and there is no net directionality of the trait mean (Hunt Reference Hunt2006, Reference Hunt2008b). The difference between stasis and an unbiased random walk is that some drift is allowable in an unbiased random walk, whereas in stasis, changes in the mean occur around a fixed point (Hunt Reference Hunt2008a). In this modeling framework, punctuated change occurs when evolutionary change is not homogenous; rates vacillate between stasis and rapid change (Hunt Reference Hunt2008b). While, undoubtedly, these modes are not mutually exclusive in the fossil record, determining the mode of morphological change of C. andinus can help to constrain the factors that may be affecting the evolution of the species. These drivers of change in the fossil lineage can be discerned by a comparison of independent records of environmental change to the fossil trajectories.
One of the fundamental features of an organism, body size, is often influenced by environmental change. Known environmental drivers of diatom cell size include instability in environmental conditions (Verdy et al. Reference Verdy, Follows and Flierl2009), energetics of lake mixing (Winder et al. Reference Winder, Reuter and Schladow2009), and nutrient availability (Litchman et al. Reference Litchman, Klausmeier and Yoshiyama2009). Our analysis focuses on the change in body size of C. andinus, and we use the likelihood framework to test which of the modes of evolution best characterize the trajectory and determine how the long-term trends in body size relate to known climate and environmental change in the Lake Titicaca region.
Background
Lake Titicaca is located on the border of Bolivia and Peru in the Andes at 3810 m above sea level (Fig. 1). The lake consists of two basins, the northern basin (Lago Grande), which is much larger (7131 km2) and deeper (maximum depth of 284 m), and the smaller (1428 km2), shallower (maximum depth of 42 m) southern basin (Lago Huiñaimarca). When lake levels are low, the two basins become separated at the Straits of Tiquina, which has a sill ~25 m below the outlet level. Dramatic lake-level lowering has occurred multiple times during the late Quaternary in response to climate change (Baker et al. Reference Baker, Seltzer, Fritz, Dunbar, Grove, Tapia, Cross, Rowe and Broda2001; D’Agostino et al. Reference D’Agostino, Seltzer, Baker, Fritz and Dunbar2002; Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007). In addition, tectonic or geomorphic change may have increased basin depth over time (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012).
Map of Lake Titicaca with depth contours of Lago Grande and Lago Huiñaimarca (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). The core material used in this study is from location 2B (star).

A 136 m drill core was recovered from 235 m water depth in the northern basin of Lake Titicaca in 2001. The recovered sedimentary record of the lake extends to ~400 Ka (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007). Global cold stages in Lake Titicaca were periods of lake-level highstands because of the positive water balance, likely resulting from decreased evaporation rates and increased precipitation, whereas the global warm interglacial intervals were registered as relative lowstands (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007). During these dry interglacial intervals (i.e., marine isotope stages [MIS] 1, 5e, 7, and 9 [Fig. 2]), Lake Titicaca became a closed basin, salinity increased, calcium carbonate was preserved on the lake floor, and saline and benthic diatom taxa dominated over freshwater planktic forms (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007, Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). Although, benthic taxa dominated the flora during interglacial periods prior to MIS 5, they were also dominant during the glacial interval MIS 8. This early exception has been attributed to changes in lake depth as a result of tectonic processes, such that the ancient lake basin was shallower than the deeper modern basin of Lake Titicaca (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012).
From left to right. The percent abundance of Chaetoceros muelleri, a diatom indicative of high saline concentrations. The relative abundance (percentage of total diatom flora) of lumped benthic species and percent calcium carbonate, which are both indicators of low lake levels. Cyclostephanos andinus sensu stricto and Cyclostephanos andinus sensu lato from the Lake Titicaca drill core (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). Images on the far right show some of the morphological variability seen in C. andinus.

During the ~400 Kyr represented by the drill-core record, the C. andinus species complex exhibits considerable morphological variation (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012) that is characterized by changes in valve size, valve undulation, and the arrangement and density of valve-face pores called areolae. The relative abundance of C. andinus also changes substantially through time. During the late Holocene (<4 Ka) and during the late Pleistocene glacial phase, MIS 4–2, 70–11 Ka, C. andinus sensu stricto (Theriot et al. Reference Theriot, Carney and Richerson1985; Tapia et al. Reference Tapia, Theriot, Fritz, Cruces and Rivera2004) was the most abundant diatom in the assemblage (Fig. 2) (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007, Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). Based on multiple lines of evidence, this species is characteristic of times when Lake Titicaca is deep, fresh, and overflowing. During the Eemian interglacial (MIS 5e; ~140 to ~115 Ka), no C. andinus morphotypes were part of the diatom flora, because lake levels were ~200 m lower than during major glacial periods, and therefore the abundance of shallow-water habitat and its proximity to the drill location favored the dominance of benthic taxa (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007, Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). From ~220 to ~140 Ka, both C. andinus sensu stricto and C. andinus sensu lato (i.e., all other morphotypes in the C. andinus species complex that were not strictly considered C. andinus by Fritz and others [Reference Fritz, Baker, Tapia, Spanbauer and Westover2012]) co-occur in the diatom flora (Fig. 2). Finally, prior to 220 Ka, C. andinus sensu lato no longer co-occurs with C. andinus sensu stricto, and the diatom assemblage was dominated by benthic taxa, again suggesting shallower lake levels. Unlike many other ancient lake diatom assemblages, Lake Titicaca maintains one planktic species or species complex throughout multiple glacial–interglacial cycles of the last ~0.5 Ma, although that group changed morphologically. In contrast, in the high-latitude Lake Baikal and Lake El'gygytgyn, for example, extinctions of taxa occur during stadial periods, followed by speciation in the interstadials (Edlund and Stoermer Reference Edlund and Stoermer2000; Khursevich et al. Reference Khursevich, Karabanov, Prokopenko, Williams, Kuzmin, Fedenya and Gvozdkov2001; Snyder et al. Reference Snyder, Cherepanova and Bryan2013 ).
Methods
Diatom Sampling
We carried out morphometric analysis on those frustule characteristics that are diagnostic of C. andinus (Tapia et al. Reference Tapia, Theriot, Fritz, Cruces and Rivera2004) from the Lake Titicaca drill-core record. This analysis included measurements of valve-face diameter (VFD), uniseriate areolae distance (UAD), biseriate areolae distance (BAD), fascicles per valve (FPV), and areolae density per 5 microns (APF) (Fig. 3). Initially, we sampled the sediment from the drill core at an interval of approximately 3000 years. This resulted in 100 samples; 55 of these samples had well-preserved measurable C. andinus valves. For samples that did not contain C. andinus, either C. andinus was not present in the flora (i.e., past lake levels were low, and the flora was dominated by benthic species), or C. andinus was not well preserved (i.e., past lake levels were very high, and valves were partially dissolved). For each sample (time step), we photographed and measured 50 well-preserved C. andinus frustules using a Zeiss Axioskop 2 plus with an AxioCam under 1000× magnification with AxioVision, Release 4.8 software. Subsequently, we carried out a higher-resolution analysis, measuring VFD every 20 cm of the 136 m core on any sample that contained well-preserved C. andinus samples. This resulted in 266 samples analyzed for C. andinus size from the ~400 Kyr record from Lake Titicaca, with an average temporal resolution between samples of ~1500 years. We used an age model established by Fritz and others (Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007, Reference Fritz, Baker, Tapia, Spanbauer and Westover2012) for the Lake Titicaca drill-core record, which was constructed using a combination of radiocarbon ages from the total organic carbon fraction of the bulk sediment, from U-series dating of discrete aragonite laminae, and by tuning peaks in the calcium carbonate sediment record to the Vostok CO2 record (Fig. 2).
Diagram of the morphometric measurements (Tapia et al. Reference Tapia, Theriot, Fritz, Cruces and Rivera2004) that were used to evaluate phyletic evolution of Cyclostephanos andinus from the drill-core record. APF, areolae per 5 microns; BAD, biseriate areolae distance; FPV, fascicles per valve; UAD, uniseriate areolae distance; VFD, valve-face diameter.

Morphometrics
To determine whether any set of morphometric measurements differentiates C. andinus morphotypes, we used principal components analysis (PCA) to condense the multiple measurements into a reduced set of components that describes the majority of the variance in the data set. We separated C. andinus into periods of time defined by MIS, because prior work showed major changes in lake level and lake chemistry (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007) and in C. andinus morphotypes (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012) in concert with these established intervals of global climate variation. We visualized the separation of MIS by using convex hulls in PCA. We also created a Pearson correlation matrix in the ‘Hmisc’ package in R to determine which of the morphometric measurements are highly correlated. In addition, we analyzed the size distribution of populations in individual samples to determine whether one or multiple populations characterize a given sample. Diatom populations commonly have a composite right-skewed unimodal cell-size distribution due to size diminution of multiple normally distributed populations through asexual reproduction (Jewson Reference Jewson1992; Jewson and Granin Reference Jewson and Granin2014). Size is reset to a maximum value through sexual reproduction after a minimum threshold is exceeded (Round et al. Reference Round, Crawford and Mann1990). In centric diatoms, this size reduction will reduce an initial valve size to 30–50% of its original valve size (Drebes Reference Drebes1977). Thus, a single population in a sample will show a single right-skewed distribution. In contrast, if the distribution is not right-skewed or there are multiple modes, this suggests that there are multiple populations present during that period. We used OriginPro 2015 to graph the size distribution from each sample that had VFD measurements (n = 266) and ran PCA in R, Version 3.2.2, on the coarse-resolution data set (six variables, n = 55).
Modeling Change in Body Size
We employed a log-likelihood approach to model the mode of body-size change of C. andinus over the history of the Lake Titicaca drill-core record (266 time steps). Here the null model is a random walk that considers two parameters—trait mean (i.e., the changes of the trait) and trait variance (i.e., the volatility of the trait) (Hunt Reference Hunt2006). We used PaleoTS, Version 0.5-1, in R, Version 3.2.2, to compare the goodness of fit of multiple modes of evolution: the general random walk (directional change), unbiased random walk (nondirectional change), stasis, and unsampled punctuated change (Hunt Reference Hunt2006, Reference Hunt2008b). We evaluated the support of each model using information criteria, specifically Akaike information criterion (AIC) and Akaike weights from the resulting log-likelihood scores (Akaike Reference Akaike1974).
Results
Cyclostephanos andinus Morphology
Many C. andinus measurements are highly correlated with size. VFD, UAD, BAD, and number of fascicles are all highly correlated and increase as overall size increases (Table 1). There are very weak correlations between metrics that vary with size (VFD, UAD, BAD, and FPV) and areolae density (APF). The first component from the PCA explains 97% of the variability in the morphological data set, the second component 2%, and the remaining three components explain less than 1% of the variation in the morphological data set. The strongest loading on PC1 is FPV, and on PC2 it is VFD. We plotted components 1 and 2 of the PCA and found that all MIS stages overlap in principal component space (Fig. 4), indicating that C. andinus morphotypes cannot be separated using the species characters quantifiable with light microscope. Preliminary scanning electron microscopy has also been uninformative in delimiting species, and we found that modern C. andinus varied very little from some of the C. andinus from the basal sediments (Supplementary Information 1). However, the PCA does reveal some patterns in the data, specifically that MIS 8 and 9 have constrained dispersion on axis 1 in principal component space, whereas, MIS 1–7 span a much larger range in morphological space (Fig. 4).
Principal components 1 and 2 from the PCA of Cyclostephanos andinus morphological traits from the low-resolution analysis of Lake Titicaca sediment samples. Convex hulls group the samples into marine isotope stages (MIS).

Correlation matrix of the morphometric measurements of Cyclostephanos andinus. The lower half of the matrix contains correlations and the upper half of the matrix reports p-values of the correlations. APF, areolae per 5 microns; BAD, biseriate areolae distance; FPV, fascicles per valve; UAD, uniseriate areolae distance; VFD, valve face diameter.

Size distributions of diatom populations are unimodal through most of the diatom record of Lake Titicaca, except during periods where C. andinus sensu stricto co-occurs with the C. andinus sensu lato (Fig. 2). After ~100 Ka and before ~200 Ka, valve-size distributions are generally right-skewed and contain one mode. At ~200 Ka, right-skewed valve-size distributions become less common (Fig. 5; Supplementary Information 2). Most size distributions that are indicative of more than a single population occurred between ~200 and ~100 Ka, which coincides with the prior delineation of multiple morphotypes of C. andinus (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007, Reference Fritz, Baker, Tapia, Spanbauer and Westover2012).
Right-skewed unimodal valve-size distribution from ~54,000 yr BP and bimodal valve-size distribution from ~180,000 yr BP.

Phyletic Evolution Modeling
The punctuated change model with two shifts in morphological trait best fit the observed data of C. andinus body size through time (Table 2). All other models—random walk, directional change, stasis, and punctuated change with one shift—had higher AIC scores, and when all of the models are weighted, punctuated change with two shifts has the highest support (Table 2). Therefore, the other models do not explain the data as well as punctuated change with two shifts. The two major shifts in size occurred at ~ 220 Ka (during MIS 7) and ~70 Ka (just after MIS 5) (Fig. 6). During the first period of stasis (~340 and ~220 Ka), mean VFD was smaller and much less variable than during the following period of stasis (~220 and ~70 Ka). The final period of stasis (~70 Ka to present) had an intermediate mean VFD in comparison with the prior two periods of stasis but retained the high variability characteristic of the period of stasis immediately preceding it.
Mean valve-face diameter of Cyclostephanos andinus for the history of the Lake Titicaca drill-core record. Marine isotope stages (MIS) are labeled above the timeline. Shaded area indicates the 95% probability interval for the best-fit punctuated change model.

Model results for the size evolution of Cyclostephanos andinus. Models tested were unbiased random walk (URW), generalized random walk (GRW), stasis, punctuated change with one shift (Punc-1), and punctuated change with two shifts (Punc-2). Highest log likelihood (logL) and lowest corrected Akaike information criterion (AICc) scores indicate the best-fitting model. K is the number of parameters in the model. Model weighting (w) shows model support.

Discussion
We found little evidence to support separating the C. andinus species complex into multiple species from the PCA of diagnostic characters apparent in light microscopy (Fig. 4). Preliminary studies of fine features of C. andinus using scanning electron microscopy also failed to provide evidence for species delimitation (Supplementary Information 1). The only statistical line of evidence for multiple populations of C. andinus is the bimodal valve-size distributions that are occasionally present from ~100 to ~200 Ka (Fig. 5, Supplementary Information 2). However, it is not evident that multiple populations indicate multiple species. Morphological work suggests that these multiple populations represent a gradient of morphological variation within the same species. Alternatively, the populations may be semicryptic species, in which size is the only discernible morphological characteristic that differentiates individuals. Regardless of whether these are multiple populations of the same species or closely related semicryptic species, they apparently were able to co-occur in Lake Titicaca during this period through spatial, temporal, or niche separation. Some type of ecological separation is required to retain multiple populations, because without separation (i.e., the driver of divergent selection and barrier to gene flow), these very closely related populations and/or semicryptic species would interbreed and fuse (Rundle and Nosil Reference Rundle and Nosil2005). The partition (spatial, temporal, or otherwise) that kept these two populations separate for ~100,000 years likely disappeared at ~100 Ka, because after this time period, the valve size of C. andinus returned to a persistent unimodal right-skewed distribution. An alternative explanation for a bimodal size distribution of C. andinus populations is preferential grazing of a particular size class (which itself could lead to selection pressure and evolution). Modern Lake Titicaca contains zooplankton that preferentially feed on specific size classes of seston (Dejoux and Iltis Reference Dejoux and Iltis1992); however, paleoecological studies of grazers have not been carried out on the drill-core record from Lake Titicaca.
While deriving process from pattern is complicated in the study of phyletic evolution, we can glean some information from model fitting. Punctuated change is the only nonhomogenous process of phyletic evolution tested here (Hunt Reference Hunt2008a). This makes punctuated change fundamentally different from all other modes of evolution, in that periods of stability are interspersed with rapid change. Nonlinearities and discontinuities are important hypotheses of how complex systems operate. In ecology, ecosystems can move abruptly from one stable regime to an alternate stable regime (Scheffer and Carpenter Reference Scheffer and Carpenter2003). Additionally, research on modeling genomic architecture has shown that gradual adaptive change acts as a driver to produce alternative stable regimes in diverging populations during the course of speciation (Flaxman et al. Reference Flaxman, Wacholder, Feder and Nosil2014). Nonlinearities appear to be an essential aspect of the functioning of complex systems, and characterizing the environmental dynamics that occur during shifts and periods of stasis can help to illuminate the process that underpins the pattern of evolution.
Modeling the body-size changes of C. andinus revealed two major shifts in valve size with intervening periods of stasis during the lake’s history (Fig. 6). The first shift in VFD occurred in MIS 7, during an interval when it has been hypothesized that Lake Titicaca’s basin evolved from a shallow basin to its modern larger and deeper mode (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). Prior to this time, C. andinus was much smaller, and the variability of VFD is much more constrained compared with the other two periods of stasis. This smaller valve size coincides with higher δ13C values, indicative of a macrophyte-dominated lake (Cross et al. Reference Cross, Baker, Seltzer, Fritz and Dunbar2000; Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007) and a period of abundant and diverse benthic diatoms in the flora, both lines of evidence that support an inference of low lake levels (Fig. 2) (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). The second shift in VFD occurred just after MIS 5, the most saline interval of the last 400 Kyr, as the result of a long period of reduced effective moisture and closed-basin conditions (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007). This is followed by the most recent period of stasis, which spans the last glacial/interglacial cycle (MIS 1–4) and contains C andinus valves with intermediate mean size with increased variability just prior to the Holocene.
The earliest jump in VFD is coincident with the transition from a relatively shallow lake to a deeper basin (Figs. 2, 6), possibly associated with regional tectonic tilting that may have raised the elevation of the outlet of the lake (Fritz et al. Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). A deeper lake may have a greater diversity of habitats available to diatom populations, dependent on the morphometry of the basin (Stone and Fritz Reference Stone and Fritz2004). This habitat variation could explain both the high variability in VFD and the bimodal distributions of size of C. andinus during the second period of stasis. This pattern is also consistent with the observation that planktic diatoms can rapidly evolve given new habitat (Theriot et al. Reference Theriot, Fritz, Whitlock and Conley2006). This deepening would provide not only new spatial niches to exploit but may have produced temporal separation by altering the mixing regime and thus the presence and duration of thermal stratification. This may allow different size classes to reproduce at different times during the year. In ancient Lake Baikal, two closely related species of Aulacoseira that have coexisted in Lake Baikal for more than 100,000 years (Khursevich et al. Reference Khursevich, Karabanov, Prokopenko, Williams, Kuzmin, Fedenya and Gvozdkov2001) exhibit niche differentiation by sexually reproducing and blooming during different seasons (Jewson et al. Reference Jewson, Granin, Zhdarnov, Gorbunova, Bondarenko and Gnatovsky2008). Two less likely alternatives are that the punctuation in VFD indicates the presence of a closely related species or that there was a hiatus in the sedimentary sequence giving the impression of punctuation. Invasions from neighboring sites have been posited as an explanation for punctuations in other fossil lineages (Van Bocxlaer et al. Reference Van Bocxlaer, Van Damme and Feibel2008). However, this seems unlikely in this case, because, as of yet, no large paleolake deposits have been identified that would contain the appropriate habitat for a source population of closely related C. andinus species. Another possibility is that a species closely related to C. andinus has come from within Lake Titicaca or evolved there during one of the punctuation events. While indeed there may be a hiatus in the Lake Titicaca drill-core record at ~120–80 Ka (Fig. 2), this hiatus does not co-occur with either of the punctuations, and there is no evidence of an abrupt change in sedimentation rate or pattern in the age-depth relationship (Fig. 2) in the interval associated with the earliest punctuation in VFD.
The second punctuation in C. andinus VFD occurred at the beginning of the glacial period MIS 4, just after the global warming of MIS 5 (Fig. 6), and is likely the result of deepening water and expanded habitat. After this shift, the average size of C. andinus was reduced, but it still retained high variability in valve size. The last interglacial (MIS 5) is marked by a near-complete removal of C. andinus populations from the diatom flora and an increase in benthic taxa, indicating low lake levels (Fig. 2) (Fritz et al. Reference Fritz, Baker, Seltzer, Ballantyne, Tapia, Cheng and Edwards2007, Reference Fritz, Baker, Tapia, Spanbauer and Westover2012). A lowered lake level would reduce lake surface area, which could decrease fetch length and reduce mixing depths, preferentially favoring smaller diatoms (Willen Reference Willen1991). During the last glacial cycle (MIS 4–2), C. andinus VFD exhibited a great deal of variability in size (Fig. 6). This may have been driven by variable nutrient input to the lake basin and/or variable turbidity within the water column associated with the frequent large-amplitude climate variations observed during this period (Fritz et al. 2010). Both nutrients and light penetration are known to influence diatom size (Willen Reference Willen1991), and both would have been variable as glaciers expanded and contracted in the surrounding cordillera and riverine discharges into the lake dramatically increased and decreased.
In conclusion, the size of C. andinus exhibited high variability over the last 400 Kyr in Lake Titicaca. While, cryptic species or semicryptic species cannot be ruled out, differentiation of populations is only easily achieved based on size distribution within an individual sample. Overall, the changes in body size of C. andinus are best explained by punctuated change driven by environmental dynamics, specifically changes in lake depth caused by either alteration of the basin morphometry or large-amplitude climate variability at a range of timescales. Habitat expansion, and possibly range expansion, is likely responsible for the initial increase in VFD of C. andinus driven by tectonic or geomorphic influences on Lake Titicaca’s basin. After the lake basin neared its modern configuration, climate variation drove variation in VFD and thus produced a pattern of prolonged stasis and rapid punctuation. In the future, genetic analysis of modern C. andinus from Lake Titicaca and other smaller, younger lakes in the Altiplano that also contain C. andinus may help to inform the history of this endemic lineage of diatoms.
Acknowledgments
We thank Maddie DeVeney, Sara Heine, Jacob Morehouse, and April Thalken for assisting in documenting Cyclostephanos andinus from the Lake Titicaca drill-core material. We also thank Jeffery Stone, Yanbin Lu, Karlyn Westover, Vicky Chraibi, and Marie Weide for insight into the development of the project, and Caroline Burberry and Scott Keely for use of their computer hardware. Comments from Ed Theriot, Craig Allen, Diana Pilson, David Watkins, Ross Secord, Gene Hunt, and one anonymous reviewer greatly improved this paper. This research was funded by the National Science Foundation Integrative Graduate Education and Research Traineeship (IGERT) program (NSF 0903469), the Sedimentary Geology & Paleobiology program (NSF 1251678), and the Frontiers in Earth Systems Dynamics program (NSF 1338694). The Geological Society of America, the Paleontological Society, the American Microscopical Society, and a University of Nebraska Presidential Graduate Fellowship also provided support for this project. This contribution is in partial fulfillment of the requirements for a doctoral degree for T.L.S. at the University of Nebraska–Lincoln.
Supplementary Material
Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.rf5jp







