1 Materials and methods 1.1 Sampling floras

For purposes of this study, our sampling units were fossil floras within the Eocene and Oligocene epochs. We defined a fossil flora as a list of macrofossil species representing one collection effort for a single well-defined stratum and geographic location (see Greenwood, 1991; DiMichele et al., 2004). Therefore, we treated collections from the same stratum at different geographic locations as representative of different fossil floras. Arguably, this resulted in some lack of independence among our fossil floras, especially among those representing geographically close locations of the same stratum. However, most studies across fossil localities must make similar decisions on what constitutes a sampling unit, and such decisions may run some unavoidable risks of either introducing unwanted autocorrelations or combining geographically distinct, botanically unrelated floras (Barghoorn, 1951; Crane and Lidgard, 1989; Greenwood, 1991).

We obtained lists of published fossil floras representing the Eocene and Oligocene of North America by using four literature sources, because no single literature source provided a comprehensive list (see Wing, 1987; also Powell, 2009 and Alroy et al., 2008 regarding comprehensiveness of the Paleobiology Database at https://paleobiodb.org/#/). Our four sources were: (1) Penhallow (1908), (2) Hollick (1936), (3) Barghoorn (1951), and (4) Wolfe et al. (1998). Each source indexed fossil floras, and for each flora, we assessed the availability of its species list (i.e., digitally or in print; from the original source or reprinted) and its consistency with our definition of a fossil flora. For the fossil floras from Penhallow (1908) and Hollick (1936), we used more recent publications to assign the fossil floras to epochs (namely MacNeil et al., 1961; Nokleberg et al., 2000) and retained only those of Eocene and Oligocene age. For all floras, we verified their geological age using Geolex (http://ngmdb.usgs.gov/Geolex/search) and a survey of the literature, and we assigned each Eocene flora to intervals comprising early, middle or early, which correspond to the Ypresian (56.0-47.9MYA), Lutetian/Bartonin (47.9-37.8MYA), and Priabonian (37.8-33.9MYA) ages, respectively. We performed our downstream analyses using the most current and widely accepted epoch and interval for each flora, but we present date ranges from Geolex or found among the literature in Appendix 1 to show where disagreements currently exist. Notably, future refinements in the geological age of the floras could affect our outcomes and necessitate updates to our analyses. In total, we recovered 201 fossil floras, of which 172 were of Eocene age and 29 represented the Oligocene (Table 1; Appendix 1).

We obtained digitized lists of genera present in the floras from the Paleobiology Database (https://paleobiodb.org/#/), by applying optical character recognition (OCR; online tool at http://www.onlineocr.net/) to.pdf files of species lists, or via manual data entry. We retained the genera as determined by the original authors, and we did not perform taxonomic reconciliation. Taxonomic reconciliation involves standardizing a dataset by applying names from a single source, or a few sources (Isaac et al., 2004). Previously, the effects of taxonomic reconciliation on genus richness were tested by Wagner et al. (2007) on a dataset comprising fossil marine mollusks. The authors obtained data records for mollusks from the Paleobiology Database and reconciled generic names by uniformly applying their own expertise and the most current literature (Wagner et al., 2007). Wagner et al. (2007) found that genus richness based on the raw and reconciled lists were similar for the three geochronological time units included in their study. Similarly, the negligible effects of reconciliation on richness were also demonstrated by comparing studies on trilobites (Foote et al., 2007). Thus, we used the generic names from the original sources to calculate genus richness of each flora.

We obtained the modern geo-coordinates of each fossil flora. Initially, we estimated the geo-coordinates for the floras using GeoLocate software (Rios and Bart, 2010) in its web application mode. GeoLocate works by parsing strings of location information, which we supplied for each flora from the original or a secondary source. In cases where GeoLocate could not resolve geo-coordinates of the fossil floras or the estimation seemed dubious, we found the coordinates using Google Maps (https://www.google.com/maps), Google Earth Pro (https://www.google.com/earth/), ME Home Town Locator (http://www.hometownlocator.com/), or American Century Farms Finder (http://www.agricultureslastingheritage.org/american-century-farms-finder/).

We used the modern geo-coordinates of the fossil floras to find their paleolatitudes. For floras east of the Rocky Mountains, we obtained paleolatitudes by using the Paleolatitude Calculator (http://www.paleolatitude.org/). The calculator accepted as input the modern coordinates, an age in millions of years, and a paleomagnetic reference frame. For age, we used the mid-point ages of the early and middle Eocene and the Oligocene, and these are, respectively, 51.9MYA, 42.8MYA, and 28.4MYA. There were no late Eocene floras east of the Rocky Mountains. Within the Paleolatitude Calculator, we applied the default paleomagnetic reference frame. For floras within and west of the Rocky Mountains, there were no available paleomagnetic reference frames in the Paleolatitude Calculator. For these, we found paleolatitudes using the Paleobiology Database, except for some floras from Alaska, for which the database did not have records. For the fossil floras from Alaska, we found paleolatitudes using printed maps for the Eocene and (shown on 50my map; Irving, 1979) and Oligocene (shown on 25my map; Smith, 1981) epochs. For a small subsample of floras occurring east of the Rocky Mountains, we compared coordinates derived from the maps, the Paleobiology Database, and the Paleolatitude Calculator and found them to be highly similar (i.e., ±1° or less; data not shown).

1.2 Data analyses

We used quadratic regressions to test for linear and curvilinear trends in the latitudinal gradients of vascular plant genera within the Eocene and Oligocene epochs as well as for the early and middle intervals of the Eocene. The late interval of the Eocene possessed too few samples (n=4) for meaningful analysis, so we performed an additional analysis in which we combined it with the middle Eocene. We applied quadratic regression, because curvilinear latitudinal gradients have been detected in some extant organisms, and we wanted to consider this possibility (Lyons and Willig, 1999; Sax, 2001; Gaston and Spicer, 2013).

We preferred not to remove floras from our dataset for having low genus richness, especially because there was no clear, nonarbitrary way to do this. In particular, we could not know the minimum number of genera required for a significant reduction in sampling biases or even if low genus richness introduced such biases in the first place. In fact, some authors tentatively propose that there the relationship between observed and actual number of genera within a fossil flora is such that floras with low observed richness should not bias the data more than floras with high observed richness (Krasilov, 1975). Nevertheless, we tested the effects of low genus richness on our regression outcomes by performing additional analyses comprising only floras having 2+ and 3+ genera. We performed these analyses for the Eocene (all intervals combined) and Oligocene datasets.

For comparison to the fossil floras, we assessed genus richness of modern North American floras north of Mexico using the Floras of North America database (http://botany.okstate.edu/floras/; also in peer reviewed articles, e.g., Palmer et al., 2006; Qian et al., 2007). The Floras of North America database comprises detailed records of published floras, which represent approximately 200 years of botanical exploration and were vetted for minimum information standards (i.e., Palmer and Richardson, 2013). From the database, we obtained genus richness and latitudinal mid-point for 3903 floras ranging in size from 3.0 × 10^-3 ha to 9.9 × 10⁸ ha latitudes of 24.6-81.8°N. We analyzed the latitudinal gradient of the modern floras using quadratic regressions as we did with the fossil floras. We also performed a multiple linear regression with latitude and area, which is expected to have a strong effect on richness but which cannot be accounted for in fossil floras. We also accounted for the effect of area on genus richness among the modern floras by graphing the residuals of the genus richness versus area relationship. This allowed us to visualize what part of remaining variation in richness could be attributed to latitude.

For all regressions, we used the log transformations of genus richness and, where applicable, of area. This transformation is consistent with the commonly-applied Arrhenius power-law model (Arrhenius, 1921, 1923), which has been shown to yield the best fit of richness data at larger geographic scales and to reduce heteroscedasticity (Hopkins, 1955; Whittaker, 1972; also at smaller scales, see Fridley et al., 2005). Hereafter, all reference to genus richness and area as variables in our regressions analyses assume that these variables have been log transformed.

2 Results

Our 201 fossil floras used for data analyses comprised 172 floras representing the Eocene epoch and 29 representing the Oligocene epoch (Table 1). The fossil floras covered a geographic extent from the southern United States to Alaska, though there were no Oligocene floras from Canada (Fig. 1; also see Berry, 1925). Within the Eocene, there were 139 floras representing the early interval, while the middle and late intervals had 29 and 4 floras, respectively (Table 1). Twenty Eocene floras and eight Oligocene floras possessed one genus each, while 19 Eocene floras and one Oligocene flora had two genera each (Appendix 1). The most diverse Eocene flora comprised 99 unique genera (Wilcox formation; Puryear, Henry County, Tennessee; Berry, 1930) and the most diverse Oligocene flora consisted of 48 unique genera (John Day formation; Fossil locality; Meyer and Manchester, 1997) (Appendix 1).

Our quadratic regression analyses of 172 Eocene and 29 Oligocene fossil floras, respectively, yielded significant quadratic coefficients and were, therefore, better fits for the data than the linear models (Fig. 2A and B). In particular, our quadratic regressions revealed hump-shaped trends in latitudinal richness for both epochs (p=0.0136, p < 0.0001 for the quadratic coefficients for the Eocene and Oligocene, respectively; Fig. 2A and B). For the Eocene, the relationship was weak (R²=0.1907), and genus richness peaked at 45.5°N (Fig. 2A). For the Oligocene, the quadratic relationship was strong (R²=0.7853), and richness peaked at 48.6°N (Fig. 2B). Our regression analyses that excluded floras with low richness showed the same significant trends with similar peaks in richness (Fig. 2C-F).

We obtained mixed results for the time intervals of the Eocene. For the early Eocene interval, we found a weak, significant curvilinear trend (R²=0.0849, p=0.0098 for the quadratic coefficient; Fig. 3A). For the middle Eocene, we found no significant linear or quadratic trends (Fig. 3B). Floras representing the middle + late Eocene interval showed a significant linear trend with increasing richness towards the poles (R²=0.3040; Fig. 3C).

Modern floras exhibited a significant quadratic relationship (p < 0.0001 for the quadratic coefficient) between genus richness and latitude and a significant linear relationship between genus richness and area. However, the relationship was weak (R²=0.0384), and richness peaked at 31.5°N (Fig. 4A). Richness was much more strongly correlated with area than with latitude for the modern floras (Fig. 4B). However, we detected latitudinal trends for the modern floras even when we performed the multiple regression with area as a term (Table 2) and when we graphed the residuals of the genus-area relationship as a function of latitude (Fig. 4C).

3 Discussion

The paleobotanical record must be applied cautiously to testing hypotheses about the latitudinal gradient of vascular plants because of the high potential for biases in the data. Among the probable sources of biases are the effects of area and time, which usually cannot be assessed for fossil floras. Area is the strongest known predictor of taxonomic richness in modern floras (Watson et al., 1835; Connor and McCoy, 1979, 2001), and this is evident from our analysis of area as a predictor of genus richness (Fig. 4B). Despite the importance of area for predicting richness, the area represented by a fossil flora usually cannot be known. Similarly, the time of accumulation for a fossil flora may not be known, but the relationship between taxonomic richness and time of accumulation is nether trivial nor linear over fine and broad time scales (Johnson, 1960; Preston, 1960; Nee et al., 1992; Rosenzweig, 1995; Fine and Ree, 2006; White et al., 2006; Gaston and Spicer, 2013). Stated another way, a fossil flora may represent several millions years of accumulation with many opportunities for the introduction or evolution of additional taxa (Behrensmeyer, 1982; Cross and Taggart, 1982). Another source of biases is sampling. One type of sampling bias arises due to taphonomy, or the preferential preservation of some taxa over others within floras (Delevoryas, 1962; Greenwood, 1991; Willig et al., 2003). Sampling bias may also result from unequal collection efforts; that is, uneven distribution of investigated floras with published taxonomic lists (Palmer et al., 2002; Zhao et al., 2006). Although we assert the importance of caution in the interpretation of our results, we also believe that similarities in latitudinal trends between paleo (Fig. 2A-B) and modern floras (Fig. 4A) may indicate that we have detected a real pattern, rather than artifactual one resulting from biases.

Nevertheless, taphonomic biases are particularly noteworthy for interpretations of richness in fossil floras and can lead to misinterpretation in several ways. One result of taphonomic bias is that most preservation of plant remains occurs within and around lake and river systems (Greenwood, 1991; Behrensmeyer, 1992). Dry systems, such as beach-front plant communities, are not well represented in our data, because they are not well represented in the paleobotanical record overall (Behrensmeyer, 1992). Therefore, our results may best represent the genus richness gradient of plant communities that occurred near water, and other systems could have had different arrangements of biodiversity during the Eocene and Oligocene of North America. Additionally, preservation potential of individual plant parts is strongly biased towards hard structures, such as hard fruits or rhizomes, and abundant parts, such as leaves of trees (Greenwood, 1991; Spicer, 1991). Plants that have few hard or abundant parts, such as herbaceous dicots and most monocots, have much more limited preservation potential (Daghlian, 1981; Spicer, 1989). Thus, while plants comprising mostly soft parts do occur in our taxonomic lists (e.g., Sparganium L. and Caulinia Willd.; Appendix 1), but they are probably underrepresented so that our findings may be most applicable to woody genera.

These preliminary results show that genus richness is highest in the mid-latitudes of North America and declines towards the border with Mexico and the Gulf Coast and towards northern Canada. This pattern is robust to floras with exceptionally low observed richness of 1-2 genera (Fig. 2). The pattern is also apparent within the early Eocene (Fig. 3A). The middle and middle + late Eocene exhibit a linear trend of increasing richness with increasing latitude (Fig. 3B and C), but the trend most likely reflects reduced sampling of high latitude floras in these time intervals rather than a real shift away from the curvilinear pattern of biodiversity arrangement.

Mean annual temperatures may explain the hump-shaped latitudinal richness gradient for Eocene fossil floras and why the peak in genus richness during the Eocene of North America was further north than for modern floras (Fig. 2A). Mean annual temperature is likely to be a strong direct or indirect (e.g., via correlation with other temperature variables) driver of the latitudinal gradient in vascular plants; warmer temperatures may facilitate greater productivity and/or speciation (Jansson, 2003; Hillebrand, 2004; Mittelbach et al., 2007; Qian, 2013; Gillman and Wright, 2014; Kerkhoff et al., 2014). Eocene temperatures were not only higher than in present day (Axelrod, 1966; Graham, 1993, 1999; Zachos et al., 2001), they also probably followed a curvilinear trend in North America, rather than a linear one (Sloan and Barron, 1992; Greenwood and Wing, 1995). In particular, Greenwood and Wing (1995) used physiognomic and taxonomic methods (CLAMP described in Wolfe, 1993 and a palm-based nearest living relative method, respectively) to estimate mean annual temperatures for many Eocene sites in North America and found a hump-shaped distribution of temperatures; namely that temperatures increased from 39°N to 43°N before declining again towards the poles (with the exception of the most southerly site, which had a very high estimated mean annual temperature). This is corroborated by a meta-study of mean annual temperatures from Eocene fossil plant localities, which showed higher temperatures in central, continental North America than along its margins, including its southern margin (Sloan and Barron, 1992). Thus, the curvilinear Eocene temperature gradient may be sufficient to explain both the stronger curvilinear trend in richness during the Eocene than in the present day as well as the higher latitude peak in richness.

Our regression results for the Eocene (Figs. 2A and 3) differed from those of Harrington (2004), who examined pollen data from the late Paleocene/early Eocene of North America and detected a linear trend of decreasing richness with increasing latitude. Harrington (2004) acknowledged that pollen often represents regional, rather than local, richness (i.e., catchment areas are large), and he accounted for this by performing his analyses using pooled pollen samples representing broad regions. Pooling samples also reduced the effects of preferential preservation of wind dispersed pollen due to its greater abundance (Crane and Lidgard, 1989; Harrington, 2004). Moreover, although Harrington (2004) studied pollen in strata that transversed the Paleocene/Eocene boundary, the effects of the boundary events (e.g., mass extinction; Sharpton et al., 1992; Schulte et al., 2010) on plant community composition was probably negligible (McIver and Basinger, 1999; Wing and Harrington, 2001; Harrington, 2008). Therefore, it is unlikely that differences between our study and Harrington's (2004) can be explained by the catchment area of pollen, biases in pollen preservation potential nor Paleocene/Eocene boundary events. Plausible explanations may include sample size, which was reduced in the Harrington (2004) study by pooling and other, unexplored differences between the pollen and macrofossil records.

In contrast to the Eocene, the hump-shaped richness gradient during the Oligocene of North America (Fig. 2B) may be attributable to evolutionary and historical factors rather than climate-related, ecological ones. Specifically, the Oligocene marked a period of intensive diversification of mesic plant species within temperate zones in response to cooler global mean annual temperatures globally and increased seasonality at mid-and high latitudes and (Donoghue et al., 2001; Zachos et al., 2001; Milne and Abbott, 2002; Donoghue and Smith, 2004; Fine and Ree, 2006; Milne, 2006; Zanazzi et al., 2007; Harris et al., 2013; Kerkhoff et al., 2014; Prothero and Berggren, 2014). Temperate species richness increased, based on studies of global pollen types and North American leaf margin features, even as tropical areas grew smaller (Wolfe and Barghoorn, 1960; Mittelbach et al., 2007). Thus, evolutionary radiation of temperate-adapted plant species may explain the peak in genus richness at mid-latitudes that we observed for Oligocene floras.

Alternatively, decreased precipitation at low latitudes may explain the hump-shaped richness gradient in the Oligocene of North America (Peterson and Abbott, 1979). Prior studies in various groups of organisms including plants have shown that waterrelated variables, such as mean annual rainfall and evapotranspiration, may be strongly, positively correlated with richness, especially at lower latitudes (Currie, 1991; Rahbek and Graves, 2001; Hawkins et al., 2003). Mean annual precipitation decreased markedly in the southern United States from~150 cm/year during the early Eocene to much lower levels of~63 cm/year at the beginning of the Oligocene (Peterson and Abbott, 1979; Woodburne et al., 2009). Thus, increasing aridity during the Oligocene at the southern extent of our study region may explain lower genus richness at lower latitudes and, consequently, the peak in diversity at midlatitudes. However, aridification during the Oligocene remains poorly understood and a topic of continued debate (Roth, 1984). Moreover, relationships between water-availability and richness appear negligible or non-existent for some organisms (Pianka, 1966).

Notably, we observed that genus richness was parabolic for both past epochs and the modern floras (Figs. 2 and 4) even though climates have changed considerably from the Eocene to present (Zachos et al., 2001). Therefore, some non-climatic factor may partially explain the latitudinal gradient. One plausible, nonclimatic factor is elevation, which usually has a hump-shaped relationship with richness so that richness is greatest at roughly mid-elevations (Sanders and Rahbek, 2012; Guo et al., 2013). The effects of elevation are likely linked with climatic factors, especially decreasing temperature with elevation, as well as with nonclimatic ones such as disturbance, productivity, and area (Sanders and Rahbek, 2012). Elevation in North America increases along a longitudinal gradient from east to west (Smith et al., 2010) since at least the beginning of the Paleocene (66 MYA) (Tweto, 1975), and our sampling of fossil floras also generally increases in latitude from east to west (Fig. 1; using modern geocoordinates). Thus, the highest latitude floras in our study may also have had the highest elevations. However, one of our sources for 14 floras (Wolfe et al., 1998; source # 4; see Appendix 1) was a study of paleoelevations in western North America, and those floras exhibit locally variable paleo-elevations from < 0.1 km to 2.9 km that do not explain richness (Supplementary File 1). Moreover, elevation is unlikely to explain the quadratic pattern in richness for modern floras, for which we have robust sampling at all latitudes of eastern and western North America (map at http://plantbio.okstate.edu/floras/). Nevertheless, future studies may generate and integrate paleo-elevational information for more floras along a latitudinal gradient as well as include more high latitude floras from eastern North America.

Another, plausible, non-climatic explanation is that the trend in genus richness varies latitudinally with available continental area. In North America north of Mexico, continental area generally increases with increasing latitude, and this has been invoked in prior studies to explain why the greatest species richness in diverse temperate organisms occurs between~30°N and 40°N instead of at the lower temperate latitudes (Sax, 2001). Similar trends in temperate richness have been detected Europe, which also exhibits increasing area with increasing latitude (Sax, 2001). Larger continental area at mid-latitudes may lead to more richness in the region, which represents the species pool for floras. However, not all studies have detected relationships between latitudinal availability of continental area and richness (Janzen, 1981). Moreover, if area plays a role in limiting richness at lower latitudes, its effects at higher latitudes, where area is greater in the Northern Hemisphere, must be mitigated by climatic or other factors to explain the lower richness in those regions.

In summary, our results provide preliminary evidence of a richness gradient in the Eocene and Oligocene of North America that is more strongly parabolic than in modern times. The humpshaped trend in richness during the Eocene may have been driven by mean annual temperature, which was probably highest at mid-latitudes in North America. During the Oligocene, peak richness at mid-latitudes may have resulted from evolutionary radiations of vascular plants within the temperate zones or from increased aridity at lower latitudes. Future studies may explicitly test these hypotheses regarding the roles of climate and diversification on the latitudinal gradient in Eocene and Oligocene of North America. Preliminarily, our results suggest that the modern latitudinal genus richness gradient of vascular plants in North America, which is much less strongly unimodal than in the past and shows peak richness at a lower latitude, may have originated after the Oligocene epoch.

Appendix A. Supplementary data

Supplementary data related to this article can be found at http://dx.doi.org/10.1016/j.pld.2016.06.002.