The Species-Area Relationship in the Late Ordovician: A Test Using Neutral Theory

The fundamental biodiversity number, θ, as proposed by Hubbell, should be positively correlated with province area. Because θ can be calculated from preserved relative abundance distributions, this correlation can be tested in the fossil record for regions with known provinces. Late Ordovician (443–458 Ma) strata of Laurentia are divided into four geochemically and biologically distinct regions that reflect provinces in the epicontinental sea. We use existing and newly obtained bed-level census data to test whether Hubbell’s θ is positively correlated with the area of these four regions, corresponding roughly to the Appalachian Basin, Cincinnati Arch, Upper Mississippi Valley, and western United States and Canada. Results indicate a positive relationship between province area and θ that suggests the influence of provincial area, among other factors, on diversity. This correlation highlights the inherent link between diversity and abundance structure at local and regional scales, such that changes at one scale will necessarily affect the other. Since diversity at these smaller spatial scales is an important component of global biodiversity, determining the nature of this relationship in the fossil record has implications for understanding how diversity is assembled globally throughout the Phanerozoic.


Introduction
Hubbell proposed a neutral theory of biodiversity and the existence of a fundamental biodiversity number θ, which can be calculated from relative abundance distributions, and, which reflects diversity and evenness within an ecological community [1].This single parameter links local and regional diversity through a positive correlation with habitat area and, in doing so, implies that knowing local diversity can aid in inferring regional diversity.Because θ can be calculated from preserved relative abundance distributions [1], it is possible to test for the expected correlation between θ and province area.
The Late Ordovician epicontinental sea of Laurentia has been divided into four geochemically and biologically distinct provinces [2,3], which reflect original oceanographic differences in water properties.Since there was likely little mixing between the water masses of these regions [4], each can be treated as a single metacommunity or a province that experiences little or no biotic interchange with adjacent metacommunities [1].By dividing the Laurentian epicontinental sea in this manner, it is possible to test for the predicted linear relationship between θ and metacommunity area.

Overview of Hubbell's Neutral Theory
Hubbell's neutral theory is a dispersal-based approach to understanding how single trophic levels of ecological communities are assembled [1].Space within a community is biologically saturated, but instead of a species' ability to occupy space being dependent on interactions and niches, it is probabilistically controlled by community size, relative abundance, and dispersal limitation [1].Species abundances change in accordance with the process of ecological drift by undergoing a random walk through time [1].Under ecological drift, the abundance of species is a zero-sum game; any change to population size of one species is balanced by a population increase or decrease in another species, such that community size is constant [1].A species that occupies more than half of the space within a community is more likely to drift toward monodominance, while a species that occupies less than half of the space in the community is more like to drift toward extinction [1].
The main assumption of the neutral theory is that all species are competitively equal, such that the per capita probabilities of birth, death, migration, and speciation are the same for all individuals, regardless of life habit or ecological role [1].This assumption is controversial among ecologists, and, although neutral models successfully predict many aspects of community dynamics, e.g., [5], many have argued that the assumption of neutrality is an inaccurate oversimplification of community ecology, e.g., [6].Although neutral theory, like all models, contains inherent simplifications of natural systems, it can be used as a test and provide valuable information even if it fails [7].
Most studies of the neutral theory have focused on modern ecological communities, but the potential for paleobiological applications, although not yet widely explored, is promising [1,7,8].Because it reduces the complexity needed to understand diversity and abundance within an ecological community, neutral theory has the potential to serve as a model that can be applied to and tested in the fossil record.Only a few studies have applied the neutral theory to the fossil record, e.g., [9,10], and further confirmation is needed that ecological communities throughout the Phanerozoic conform to patterns predicted by an assumption of neutrality.This could have significant implications for a wide array of paleobiological questions [1,7,8].

The Parameters θ and m
An important aspect of neutral theory is the fundamental biodiversity number θ, a dimensionless number that, along with the dispersal coefficient (m) governs the shape of a community's relative abundance distribution [1].As a reflection of species richness (diversity) and the similarity of species' relative abundances (evenness), θ controls the length and slope of the relative abundance distribution of a community [1] (Figure 1).An increase in θ results in an increase in diversity and a decrease in monodominance as space in the community becomes occupied by more species with similar abundances [1].The relative abundance distribution increases in length and decreases in slope until θ approaches its theoretical limit of infinity, at which point, the shape of the distribution approaches a horizontal line, and all species within the community have an abundance of one [1].

Figure 1.
The effect of the fundamental biodiversity number, θ, on relative abundance distributions.Communities with a larger θ are more diverse and have greater evenness.Modified from Hubbell [1].
The dispersal coefficient (m) is the probability that a death in the local community will be replaced by an immigrant from outside the local community but within the metacommunity [1].A metacommunity is a closed system; species can migrate between local communities within a metacommunity but not into or out of the metacommunity itself.As a result, m is a measure of dispersal limitation within the metacommunity.It also indicates patchiness within a metacommunity and reflects the amount of similarity between local and metacommunity abundance structures.A local community with an m close to one experiences less dispersal limitation and will have a relative abundance distribution that is more similar to the metacommunity than a local community with an m close to zero [1].Since most local communities experience some degree of dispersal limitation, and Diversity 2013, 5 243 because species that are rare in the metacommunity have a lower probability of dispersing to a local community, the relative abundance distribution of a local community will be shorter and steeper than that of the metacommunity [1] (Figure 2).

Figure 2.
The effect of the dispersal coefficient, m, on the relative abundance distribution of a local community compared to that of the metacommunity.As m decreases, local communities (thin lines) differ increasingly from the metacommunity (bold line).Modified from Hubbell [1].

The Relationship between θ and Area
Hubbell derived θ from modeling the random walk of a dispersal-assembled community with a metacommunity size (J m ) and per-individual speciation rate (ν) [1].From these two variables, Hubbell identified a constant parameter, 2J m ν, and called it θ [1].Although this constant mathematically describes θ, it cannot be used to easily calculate θ because J m is a very large number, ν is a very small probability, and both are nearly impossible to determine in practical applications [1].In practice, likelihood methods are used to fit θ to an observed relative abundance distribution of species within a community.
Since neutral theory assumes that the space a community occupies is biologically saturated, the size of the metacommunity (J m ) is equal to the total area occupied by the metacommunity (A m ), and the mean density of individuals in that area (ρ) [1]: This equation links the biodiversity of a metacommunity to the area occupied by the metacommunity, equivalent to a biogeographic province [1].This linear relationship between diversity and province area is testable because θ can be estimated from census data, and the area of provinces can be measured from maps.

The Late Ordovician of Laurentia
In the Late Ordovician (443-458 Ma), Laurentia, the core of the continental landmass that is now North America, straddled the equator and lied entirely within tropical latitudes [11].A shallow tropical epicontinental sea, similar to the modern day carbonate platform in the Bahamas, covered most of what is now the central United States and Canada.The southeastern margin of the continent was the site of the Taconic orogeny, which created a deeper water foreland basin in the area near the modern Appalachian Mountains [12].
Bathymetric differences across the Laurentian sea coupled with its vastness created oceanographically distinct geochemical regions [2][3][4].Since these regions correspond to known faunal provinces, e.g., [13][14][15][16], they can be treated as separate metacommunities, or biogeographic provinces: the Appalachian Basin, Southern, Midcontinent, and Western provinces (Figure 3).[2] and [3], paleogeographic reconstruction from [17]).This study considers the area of only the deep subtidal facies, which for the Western Province is much smaller than total province area.The estimated extent of deep subtidal facies in the Western province is indicated by cross-hatching.Points represent sampling locations for the six data sets collected (BH: Bighorn Mountains, M: Manitoba, UMV: Upper Mississippi Valley, CA: Cincinnati Arch, AB: Appalachian Basin).

Geochemical Provinces
The Appalachian Basin, Southern, and Midcontinent provinces were defined primarily on differences in the ε Nd isotopic composition of conodonts (small phosphatic teeth of primitive fish) and the δ 13 C isotopic composition of limestone [2].The ε Nd isotopic signature of conodonts reflects the crustal age of the continental landmass that serves as the primary source of sediment to a water mass, and it can vary greatly between regions [18].The δ 13 C isotopic ratio of limestone varies only slightly in response to variations in biological productivity and burial, but can vary across carbonate platforms, e.g., [19].Holmden et al. [2] used differences in ε Nd and δ 13 C in the Late Ordovician to identify strata that record deposition within the Appalachian Basin, Southern, and Midcontinent regions.In addition, differences in δ 13 C isotope values among these regions has been linked to different rates of carbon cycling [4], further supporting the evidence for limited mixing between water masses and linking these geochemical distinctions to biological activity.That these regions had distinct isotopic compositions suggests the presence of barriers to biological migration, and thus, that these areas would have been distinct biogeographic provinces.

Biogeographic Studies Support These Regions as Provinces
Although no comprehensive studies of provinciality in the Late Ordovician seas of Laurentia have yet been conducted, studies on particular taxa indicate that the geochemically defined water masses of Holmden et al. [2] are reflected in faunal distributions.For example, the faunal distinctiveness of the Appalachian Basin Province has been recognized previously in multiple studies, including those of conodont [13,14,20], graptolite [21][22][23][24], ostracod [25], and shelly fossil [26] assemblages.Studies of solitary rugosan corals indicate a province stretching from northern Illinois southward to Missouri and southwestward to Oklahoma [27], corresponding to the Midcontinent Province of Holmden et al. [2].Rugosan corals [28] and bryozoans [16] support the presence of a province encompassing the Nashville Dome of Tennessee and the Cincinnati Arch of Ohio, Indiana, and Kentucky, similar to the Southern Province of Holmden et al. [2].Finally, the vast region of the western United States and Canada has long been recognized as a faunally cohesive region, based on bryozoans [16], nautiloids [29], rugosans [15], brachiopods [30][31][32], and a distinctive Thalassinoides ichnofacies [3,33].
The boundaries of these faunally defined provinces do not always coincide precisely with one another or with the geochemically defined water masses.Many of these discrepancies can be attributed to how provincial boundaries are drawn between sampling locations.Other differences in boundary placement likely reflect that boundaries between provinces may have been permeable to varying degrees, with some taxa more able extend to nearby regions than others.Any uncertainty in the geographic extent of provinces caused by these factors would weaken the species-area relationships examined here.

Methods
Faunal censuses were used to produce estimates of θ from each of the provinces, whose areas were measured to test for a positive relationship between θ and area.

Data Collection
To minimize variation over time, faunal samples were obtained from the C2 and C5 depositional sequences [34] (Figure 4).To minimize variation in faunal composition caused by water depth, e.g., [35,36], all samples were collected from the deep subtidal depositional environment, that is, between fair weather and storm wave base.
Six data sets were used in this study (Figure 4): three from existing data (C2 Cincinnati Arch, C5 Cincinnati Arch, Bighorn Mountains) and three from newly collected data (Appalachian Basin, Upper Mississippi Valley, Manitoba).The C2 and C5 Cincinnati Arch data sets were collected from the Fairview Formation (C2) and Liberty Formation (C5) in the Cincinnati, Ohio area [36].The Bighorn Mountains data set was obtained from the Rock Creek Beds of the Bighorn Dolomite (C5) in Wyoming [37].The Appalachian Basin data set was collected from the Reedsville Formation (C2) near Catawba and Newport, Virginia, and from Germany Valley, West Virginia.The Upper Mississippi Valley data set was assembled from the Elgin Member of the Maquoketa Formation (C5) in Fillmore County, Minnesota.The Manitoba data set was obtained from the Gunn Member of the Stony Mountain Formation (C5) just north of Winnipeg, Manitoba.Each regional data set contains multiple local samples, with new data sets containing 25-30 samples and older data sets having 40-100 samples.Locality descriptions for the new data sets are included in Appendix 1.
In most data sets, sampling occurred over multiple beds, with each faunal sample representing a single bed.Samples for the Manitoba data set, however, were obtained as individual rock slabs from a spoil pile in a quarry, and it is unclear to what extent these slabs represent multiple beds.Censuses of marine invertebrates were conducted by identifying and counting all individuals in each sample.Individuals were identified to species where possible.All taxa were counted using whole-faunal counts with a minimum number of individuals approach [35].This method calculates species abundance for brachiopods and bivalves as the total number of articulated valves, plus the greater of right or left disarticulated valves, plus one-half the number of indeterminate valves [35].Bryozoans were identified to genus when possible, but mostly by morphological form and counted in 1 cm lengths to standardize for biovolume [35].Colonial corals were also counted in 1 cm lengths, but were identified to genus.Solitary rugose corals and gastropods were identified to genus and counted as a single individual.Trilobites were counted only when identifiable, such as for pygidia and cranidia.Crinoids were noted if present but not counted since it is not possible to generate a number of individuals from columnals.
Prior to analysis, all data sets were culled to include only low-tier primary consumers to meet the neutral theory's requirement of a single trophic level [1] (Table 1), with trophic level determined from the Paleobiology Database [38].Primary consumers were selected for analysis because of their high abundance and recognizability in Ordovician samples.Focusing on only low-tier organisms, that is, those that live close to the sediment water interface, minimized any difference in nutrient access that height in the water column could provide to an organism.Samples with fewer than 24 individuals were removed from the data sets to ensure that each sample accurately represented the overall community and to increase statistical power, e.g., [35].

Relative Abundance Distributions
Two types of relative abundance distributions were calculated in R [39].A distribution that reflects the abundance distribution for local communities was calculated from the mean abundance of each species across all samples within a data set.A distribution that reflects the abundance distribution of the metacommunity was calculated from the total abundance of each species in the whole data set.Comparing the local and metacommunity relative abundances in this manner illustrates how well diversity in an average sample reflects the diversity in the province [1].

Estimation of θ
Hubbell's θ was estimated using Etienne's 2007 [40] and 2009 [41] likelihood methods.These methods calculate θ for the metacommunity and m for the local communities by using a numerical optimization over a range of possible values to determine the combination of θ and m with the maximum likelihood.These likelihood methods were performed in PARI/GP [42].
The major difference between the two methods is that the 2007 method assumes that all local communities experience the same amount of dispersal limitation, and it estimates a single m value for all local communities within the metacommunity [40].The 2009 method allows each local community to experience a different amount of dispersal limitation, and it produces a different estimate of m for each local community [41].Consequently, the 2007 method searches over a 2-dimensional parameter space, runs in a few minutes, and produces a visualizable likelihood surface.The 2009 method searches over an n+1 dimensional space where n is the number of samples, runs for a few hours for smaller data sets to several months for larger data sets, and does not result in a visualizable likelihood surface.

Estimation of Area
Province area was calculated using ArcGIS [43].The calculate geometries tool was used to measure the area of polygons of the estimated size and shape of each province on an equal-area map of North America.The approximate geographic location and extent of the provinces were determined using published maps and data from previous studies [2,3] and personal knowledge of Ordovician faunal distributions.
Although the deep subtidal environment was broadly distributed across large portions of the Midcontinent, Southern, and Appalachian Basin provinces, it is confined to a relatively small portion of the Western province.Well log data [44][45][46][47] were used to identify known subsurface locations of the deep subtidal Gunn Member and correlative strata, such as the Rock Creek beds of Wyoming [44].Additional studies, e.g., [48] of locations at which the deep subtidal are absent and of the location of anhydrites associated with the center of the Williston Basin, e.g., [47] were used to constrain the estimated extent of the deep subtidal facies to the north, south, and east (Figure 3).Because of widespread erosion of Ordovician rocks during the Devonian [46], the western extent of the deep subtidal facies could not be constrained more than the estimated western boundary of the province.

Relative Abundance Distributions
All data sets used in this study are dominated by a few abundant species (Figure 5).The metacommunity abundance distribution for the Appalachian Basin data set is the steepest and has the fewest species, and the distribution for the Bighorn Mountains has the lowest slope and the greatest species richness.In general, data sets with a larger total number of species in the metacommunity (e.g., Cincinnati Arch and Bighorn Mountains) display a long tail of rare species, which mitigates the abundance value of the most dominant species (Table 2).
The long tail of rare species likely reflects the large number of samples because increased sampling results in the -unveiling‖ of more singleton species [49].However, as long as the data set contains local samples that when combined are representative of the metacommunity [1], the effect of sample size on the estimate of  should be minimal.In all data sets except the Upper Mississippi Valley, the means of local community abundances are a good approximation of the metacommunity abundances for taxa that are neither the rarest nor the most abundant.For the Upper Mississippi Valley data set, the local community abundance distribution poorly approximates the metacommunity abundance distribution, and most local samples are much less diverse than the metacommunity.For all data sets, the most abundant taxa in the metacommunity tend to be slightly overrepresented in the local communities, whereas rare taxa and singletons in the metacommunity are often absent from individual local collections.Consequently, total species richness and evenness in local communities is less than that of the metacommunity (Table 2).

Estimation of θ and Area
Estimates of θ (Table 3) fall within a range consistent with values estimated by Hubbell for modern benthic marine invertebrate communities [1].Previously reported values of  range from 3 for a low diversity boreal community to 20 for a highly diverse tropical community.Values of θ estimated in this study are reasonable given this range, but are somewhat low given that the epicontinental sea of Laurentia was positioned entirely within tropical latitudes.
The Appalachian Basin data set has both the lowest θ and area, and the Bighorn Mountains data set has the largest θ and area (Table 3).Manitoba has an unusually low θ given the large area of the Western Province.In addition, although both the Manitoba and Bighorn Mountains data sets represent the Western province, their estimated values of θ differ.Although θ was calculated using both Etienne methods [40,41] only the results from the 2007 method [40] are used for subsequent analysis.The 2009 method [41] successfully calculated a θ for only three of the data sets (Appalachian Basin, Upper Mississippi Valley, and Bighorn Mountains).For the other three data sets, the likelihood estimation either did not finish running after three months because the data sets were too large (C2 Cincinnati and C5 Cincinnati) or produced an error (Manitoba).For the three data sets for which the multi-dimensional 2009 method [41] was successful, the estimated θ is similar to that of the 2-dimensional 2007 method [40], suggesting that the results of the 2007 method are reliable.

Likelihood Surfaces
Likelihood surfaces illustrate the uncertainty in the estimates of θ and m (Figure 6).The confidence region is smallest for the three existing data sets (C2 Cincinnati, C5 Cincinnati, Bighorn Mountains), reflecting the greater number of samples in these data sets.For all data sets, m is tightly constrained, whereas the estimate of θ has more uncertainty.The estimated values of m are less than 0.01 for all data sets except Manitoba, which has a slightly higher m (0.064).These values indicate strong dispersal limitation and little to no migration from the metacommunity into local communities [1].Such large dispersal limitation is unexpected for marine organisms, which disperse by releasing larvae into the water column.

Correlation between θ and Area
Province area and θ are weakly and positively correlated, although below statistical significance (R 2 = 0.03, p = 0.34).The trend line of the correlation between θ and area is positive (Figure 7) with a Diversity 2013, 5 252 slope of 2ρν.If the Manitoba data set is treated as an outlier and removed from analysis (see below), the correlation strengthens and becomes statistically significant (R 2 = 0.85, p = 0.02).

The Correlation between θ and Area
The weakly positive correlation between θ and area indicates that province area exerts some control on regional biodiversity.When the Manitoba data set is excluded, the correlation is substantially stronger suggesting that province area is a significant driver of regional biodiversity.
The collecting conditions in Manitoba differ from those of all other locations, which might explain the unusual estimates of θ and m for the Manitoba data set and would support treating it as an outlier.At all other sampling locations, it was possible to collect samples that were clearly from different beds.Due to limited exposure in Manitoba, samples were available only from a spoil pile in a quarry.Because of this, these samples may represent many fewer beds than the total number of samples would suggest.If all of the Manitoba samples were collected from the same bed, they would represent repeated sampling of a single local community as opposed to multiple local samples of the full metacommunity.Not sampling the full metacommunity would result in an artificially larger m and smaller θ.The Manitoba data set does have a lower θ than expected and it is the only data set to have an m greater than 0.01.
The slope of the trend line reflects a constant relationship between the density of organisms per unit area within the habitat (ρ) and the per-individual probability of speciation (ν), but it is likely that these factors vary among the metacommunities substantially enough to affect biodiversity.A stronger than expected influence of one or both of these variables would introduce scatter in the θ-area correlation and account for some of the deviation from the expected linear relationship.A large ρ or ν would result in a larger θ for a data set and cause it to plot above the trend line, while a small ρ or ν would result in a smaller θ, causing the data set to plot below the trend line (Figure 8).Differences in habitat density could easily be caused by differences in food supply or substrate.For example, a region with high productivity and more suspended organic matter could support a higher density of suspension feeders than one with lower productivity [50].This could explain why the Bighorn Mountains and Cincinnati C5 data sets, which both contain a tropical fauna, plot above the trend line while the Appalachian Basin and Cincinnati C2 data sets, which contain a cooler water fauna, plot below the line [51].
A latitudinal temperature gradient between the cooler water Appalachian Basin fauna and the equatorial Western fauna could cause variation in speciation rate, with higher speciation rates occurring in the warmer water regions, e.g., [52][53][54].However, since Hubbell's ν is a per-individual speciation rate as opposed to a per-species speciation rate as is generally measured, it is difficult to test whether this gradient is indeed a source of variation in .
An additional explanation for scatter in the data is underestimation or overestimation of province area.Overestimates of area shift points to the right of the trend line and underestimates of area shift points to the left (Figure 8).Although lithological, faunal, and geochemical data were all used to provide evidence for province boundaries, there is still a degree of uncertainty in province area mainly because burial of strata in the subsurface and erosion of strata limit the accessibility of Ordovician rocks.This could artificially weaken or strengthen the correlation between θ and area and be especially problematic for provinces in which there is a large difference between the area of the deep subtidal facies and the total area of the province.
Because sampling was limited to deep subtidal strata, uncertainty in the extent to which this habitat spanned the province could add to the uncertainty in area.For the Western province, the deep subtidal habitat is demonstrably limited to only a portion of the province, and it is this smaller area that was used in the analysis.The deep subtidal facies is likely also limited, but to a lesser extent, for the Midcontinent province, but is much more difficult to estimate.An overestimation of the deep subtidal habitat area could explain why the Midcontinent province plots to the right of the trend line.It is unlikely that possible error in the area estimate has a significant effect on the Southern and Appalachian Basin provinces because both regions have widespread and well-constrained regions of deep subtidal strata.
Additionally, the observed scatter in the relationship between  and province area could reflect non-neutral community dynamics.If interspecific competition has a significant effect on the abundance structure of the community, the assumption of neutrality would be invalid.Competition within a community could increase under environmental conditions that limit resource availability.

The Species-Area Relationship and θ
Hubbell's θ is significant for diversity studies because it links diversity across spatial scales.From the local to the global scale, the species-area relationship is triphasic, e.g, [55][56][57][58].Hubbell [1] argued that θ and m are important controls on the shape of the triphasic curve at each spatial scale (Figure 9).At the local scale, species richness increases curvilinearly with area, reflecting a sampling curve in which new species are initially encountered rapidly but then more slowly as area increases [1].At the regional scale, species richness increases log-linearly as the addition of new local communities increases the potential to add new species with each increase in area [1].At the global scale, diversity initially increases exponentially but then the slope approaches unity as metacommunities with entirely different species compositions are encountered [1].
At the local scale, θ and m have a combined effect on the species-area relationship [1].The slope of the initial part of the curve steepens with θ, such that a larger θ corresponds to a more rapid rise in the number of species as a local community is sampled.Total diversity of the local community, or alpha diversity, which occurs at the point on the curve at which the relationship switches from local to regional dynamics, is affected by both θ and m.Larger values of θ and m result in a greater alpha diversity, while lower values of θ and m will lower alpha diversity.
Diversity 2013, 5 255 At the regional scale, θ and m have inverse effects on the species area relationship [1].A larger θ will increase both the slope of the log-linear relationship (beta diversity) and the province diversity (gamma diversity).A larger m has the opposite effect, decreasing beta and gamma diversity.Total gamma diversity also corresponds to the point at which the species-area relationship shifts from regional to global dynamics, called the correlation length.This point is affected by m and ν, such that a larger ν causes a smaller correlation length and a larger m produces a larger correlation length.At the global scale, the effects of θ and m on the species-area relationship are uncertain.At this scale, species richness and area increase as multiple provinces are combined to accumulate global diversity [1].The amount of diversity added with each province depends on the degree to which provinces differ from each other.Total global diversity increases with the number of provinces and their dissimilarity from each other.
Since θ and m drive alpha, beta, and gamma diversity, diversity at local and provincial scales are linked, such that processes affecting diversity and abundance in the local community have inherent implications for diversity and abundance in the metacommunity.This link between local communities and the metacommunity reflects a continuous landscape over which local-scale diversity patterns increasingly resemble that of the metacommunity at larger spatial scales [59].
The manner in which diversity increases between these two spatial scales is important to describing how diversity is assembled globally and may help to understand the origin of changes in Phanerozoic diversity [60][61][62][63].Since the θ-area relationship exists in the Late Ordovician, and likely throughout the fossil record (though future work will need to test this), it should be possible to predict species richness of a given province from the area and habitat density of a province and an estimate of θ, as suggested by Hubbell [1]: Province area can be estimated from paleogeographic maps and θ can be estimated from a collection of local samples, so, given a method to estimate habitat density, this equation could be used to estimate species richness for any province in the Phanerozoic.

Conclusions
(1) The fundamental biodiversity number (θ) and province area are positively correlated in Late Ordovician marine communities of Laurentia, as predicted by Hubbell's [1] neutral theory.The strength of the correlation is likely affected by uncertainties in the estimates of province area and variations in the habitat density per unit area (ρ) and the per-individual speciation rate (ν).
(2) Hubbell's θ and m link local and province-scale diversity across the triphasic species-area relationship [1,59].Because of these linkages, local community diversity has implications for the diversity of provinces.This linkage could be exploited to understand causes for changes in Phanerozoic diversity.

Figure 3 .
Figure 3. Late Ordovician paleogeography of Laurentia, showing the shallow seaway (blue), land (brown), and the location of four biogeographic provinces (provinces adapted from both[2] and[3], paleogeographic reconstruction from[17]).This study considers the area of only the deep subtidal facies, which for the Western Province is much smaller than total province area.The estimated extent of deep subtidal facies in the Western province is indicated by cross-hatching.Points represent sampling locations for the six data sets collected (BH: Bighorn Mountains, M: Manitoba, UMV: Upper Mississippi Valley, CA: Cincinnati Arch, AB: Appalachian Basin).

Figure 4 .
Figure 4. Strata and time intervals sampled.Six data sets (Bighorn Mountains, Manitoba, Upper Mississippi Valley, Cincinnati Arch, Appalachian Basin) were collected from the C2 and C5 depositional sequences [34] across all four provinces.

Figure 5 .
Figure 5. Relative abundance distributions for each data set, showing the comparison between total metacommunity abundances (dark brown line) and means and standard deviations of local community abundances (tan line).

Figure 6 .
Figure 6.Likelihood contour surfaces for each data set.Maximum likelihood estimates from the 2007 method [40] are indicated with black dots.

Figure 7 .
Figure 7. Correlation between θ and province area.The dotted lines are trend lines and have a slope of 2ρν.The grey dotted line includes all data sets, and the blue dotted line excludes the Manitoba data set.

Figure 8 .
Figure 8. Causes of scatter around the trend line between θ and province area.

Figure 9 .
Figure 9.The triphasic relationship between species and area.The relationship between species richness and area varies across spatial scales and is controlled by θ and m.The effects of θ and m on the shape of the curve and local (alpha), regional (beta), and global (gamma) diversity are indicated by arrows.Based on Hubbell [1].

Table 1 .
List of taxa included and excluded in the study.

Table 2 .
Comparison between metacommunity and local community species richness, including sample sizes for each data set.

Table 3 .
Results from the estimation of θ and area.The multi-dimensional 2009 method was not able to calculate θ for some data sets due to sample size.