Quantitative Biofacies Analysis to Identify Relationships and Reﬁne Controls on Paleosol Development, Prince Creek Formation, North Slope Alaska, USA

: Late Cretaceous coastal plain deposits of the Prince Creek Formation (PCF) offer a rare glimpse into an ancient, high-latitude, arctic greenhouse ecosystem for which there is no modern analog. Here, we employ quantitative biofacies analysis to explore the spatio-temporal variability in PCF palynomorph and microbiota assemblages from nine paleosol horizons exposed along the Colville River, North Slope, Alaska. Biofacies results provide insight into paleoenvironmental controls on the coastal plain ecosystem. Cluster and ordination analyses recognize ﬁve biofacies and the following two assemblage types: (1) fern and moss dominated assemblages and (2) algae dominated assemblages. Ordination arrays biofacies along environmental gradients related to soil moisture and marine inﬂuence. Fern and moss dominated biofacies from regularly water-logged paleosols along lake and swamp margins on the lower delta plain clearly segregated from algae dominated assemblages of periodically drier levee-overbank paleosols. These results support previous interpretations from the sedimentology, paleopedology, and geochemistry of PCF paleosols that suggest that ﬂuctuations in the water table, related to seasonal river discharge and variations in topography and drainage, controlled soil development and vegetation growth across the coastal plain. This quantitative biofacies-based approach provides an independent predictive tool and cross-check for interpreting environmental conditions along any ancient coastal ecosystem. using a multivariate analytical workﬂow and biotic variation is quantiﬁed within and among PCF localities, soil horizons, and depositional environments. Environmental drivers of biofacies variability are interpreted through the integration of biotic data with published observations on paleosol sedimentology, paleopedology, and geochemistry. We suggest that this integrated approach to ecosystem

The stratigraphic and paleoenvironmental framework of the PCF coastal plain has been established through studies of continental through shallow marine deposits exposed along the Colville River [13,[21][22][23]. PCF sediments were shed off the Brooks Range to the south and west of the study area ( Figure 2) and are composed of conglomerate, sandstone, siltstone, mudstone, carbonaceous shale, coal, bentonite, and tuff. The PCF comprises the most proximal deposits of a Late Cretaceous to Paleocene progradational succession [24][25][26]  The stratigraphic and paleoenvironmental framework of the PCF coastal plain has been established through studies of continental through shallow marine deposits exposed along the Colville River [13,[21][22][23]. PCF sediments were shed off the Brooks Range to the south and west of the study area ( Figure 2) and are composed of conglomerate, sandstone, siltstone, mudstone, carbonaceous shale, coal, bentonite, and tuff. The PCF comprises the most proximal deposits of a Late Cretaceous to Paleocene progradational succession [24][25][26] and becomes younger to the north in the study area (down river) based upon the relative dating of paleosols, regional 40 Ar/ 39 AR dating [27], regional dip, palynology [28][29][30][31][32], and the presence of Paleocene (Danian) ostracods and mollusks in overlying strata [23,33,34]. Recent studies on the exposed paleosols of the PCF combined observations on their micromorphology, paleontology, geochemistry, and stratigraphic context to suggest that the PCF coastal plain represents a low gradient, muddy, ash-covered, tidally modi-  [23][24][25][26].
The purpose of this study is to quantitatively analyze palynomorph and microbiotic assemblages contained within the PCF paleosol horizons described previously by Flaig et al. [22]. Biofacies are defined using a multivariate analytical workflow and biotic variation is quantified within and among PCF localities, soil horizons, and depositional environments. Environmental drivers of biofacies variability are interpreted through the integration of biotic data with published observations on paleosol sedimentology, paleopedology, and geochemistry. We suggest that this integrated approach to ecosystem analysis improves the confidence of interpretations, suggests refinements, and aids in identifying underlying physical controls on ecosystems. These quantitative techniques can be applied broadly to examine variability among a host of geologic attributes that may not have been previously considered as candidates for multivariate analysis. The purpose of this study is to quantitatively analyze palynomorph and microbiotic assemblages contained within the PCF paleosol horizons described previously by Flaig et al. [22]. Biofacies are defined using a multivariate analytical workflow and biotic variation is quantified within and among PCF localities, soil horizons, and depositional environments. Environmental drivers of biofacies variability are interpreted through the integration of biotic data with published observations on paleosol sedimentology, paleopedology, and geochemistry. We suggest that this integrated approach to ecosystem analysis improves the confidence of interpretations, suggests refinements, and aids in identifying underlying physical controls on ecosystems. These quantitative techniques can be applied broadly to examine variability among a host of geologic attributes that may not have been previously considered as candidates for multivariate analysis.

Materials and Methods
Palynomorph and microbiota abundance data were collected from nine paleosol suc-

Materials and Methods
Palynomorph and microbiota abundance data were collected from nine paleosol successions by P.  fauna. Bulk and in situ hand samples were gathered from deeply trenched outcrops at 15 to 30 cm intervals from representative paleosols. Raco brand metal electrical boxes, measuring roughly 11.5 cubic inches in dimension, were hammered into the outcrops to collect in situ samples. Prior to use, each box's screw holder tabs were removed to make sediment collection more efficient. The boxes were oriented vertically with respect to bedding and the tops and bottoms of specimens were identified. Samples were dried out for several months and later impregnated with epoxy resin. Petrographic thin sections were later prepared by Spectrum Petrographics (Vancouver, WA, USA). Additional sample material was sent to Sue Matthews at Palynological Laboratory Services where all fossiliferous material was identified and tabulated from selected soil horizons.
exposures along the Colville River and include the Liscomb Bonebed (LBB), Kikiakrorak River Mouth (KRM), Sentinel Hill (SH), and North Kikak-Tegoseak localities ( Figure 2). No single locality exposes all nine of the paleosol horizons analyzed in this study. Figure  3 displays measured stratigraphic sections at each locality, the locations of the sampled paleosols within those successions, and their relative stratigraphic relationships to each other based on location along the Colville River, bedding dip, and dating of ash beds and fauna. Bulk and in situ hand samples were gathered from deeply trenched outcrops at 15 to 30 cm intervals from representative paleosols. Raco brand metal electrical boxes, measuring roughly 11.5 cubic inches in dimension, were hammered into the outcrops to collect in situ samples. Prior to use, each box's screw holder tabs were removed to make sediment collection more efficient. The boxes were oriented vertically with respect to bedding and the tops and bottoms of specimens were identified. Samples were dried out for several months and later impregnated with epoxy resin. Petrographic thin sections were later prepared by Spectrum Petrographics (Vancouver, WA, USA). Additional sample material was sent to Sue Matthews at Palynological Laboratory Services where all fossiliferous material was identified and tabulated from selected soil horizons.   [22] and the stratigraphic location of Paleosol 1 through Paleosol 9. Inset shows the relative stratigraphic (age) relationships of paleosols to one another as determined from regional 40Ar/39AR dating [27], structural dip, and biochronology [23,[28][29][30][31][32][33][34]. Paleosols are ordered in the inset from oldest (P1) to youngest (P9) and generally young toward the north; refer to text and Figure 1 for further details. Revised after [22].
The final dataset examined in this study includes 29 samples and 6620 fossil elements comprising 107 taxa (see Table S1 in Supplementary Material). Fossil counts vary among samples from 185 to 262 fossil elements (mean = 228 elements). Table 1 summarizes the final sample inventory parsed by paleosol horizon and environment. The sampling of depositional environments from some paleosols is incomplete because these environments were not exposed at outcrop. Fossils include in situ and reworked dinocysts and acritarchs, brackish and freshwater algae, exotic projectate pollen, lowland tree and herb pollen, bisaccate pollen, fern and moss spores, and fungal hyphae. The analyses were all carried out at the genus level. Prior to performing the statistical analyses, taxon  [22] and the stratigraphic location of Paleosol 1 through Paleosol 9. Inset shows the relative stratigraphic (age) relationships of paleosols to one another as determined from regional 40Ar/39AR dating [27], structural dip, and biochronology [23,[28][29][30][31][32][33][34]. Paleosols are ordered in the inset from oldest (P1) to youngest (P9) and generally young toward the north; refer to text and Figure 1 for further details. Revised after [22].
The final dataset examined in this study includes 29 samples and 6620 fossil elements comprising 107 taxa (see Table S1 in Supplementary Material). Fossil counts vary among samples from 185 to 262 fossil elements (mean = 228 elements). Table 1 summarizes the final sample inventory parsed by paleosol horizon and environment. The sampling of depositional environments from some paleosols is incomplete because these environments were not exposed at outcrop. Fossils include in situ and reworked dinocysts and acritarchs, brackish and freshwater algae, exotic projectate pollen, lowland tree and herb pollen, bisaccate pollen, fern and moss spores, and fungal hyphae. The analyses were all carried out at the genus level. Prior to performing the statistical analyses, taxon abundances were log transformed to emphasize the relative contributions of all taxa, effectively increasing the weight of rare elements [50]. Then, sample counts were percent transformed to alleviate any potential impact of sample size differences. Finally, Euclidean distances [50] were computed between each pair of samples and captured in a distance matrix to assess differences in their biotic composition. Two multivariate statistical methods were used to summarize and graphically display distance relationships among samples and taxa, and to interpret gradients of biotic change among PCF localities, horizons, and depositional environments; these include: (1) hierarchical agglomerative cluster analysis (HCA) with ward's method and (2) detrended correspondence analysis (DCA) ordination. These methods excel at summarizing large quantitative datasets, extracting dominant patterns, and plotting trends in clear intuitive displays. HCA and DCA have been used successfully for decades in many disciplines, but especially so in paleobiology and ecology where they reveal similarities and differences in the composition of biotas through space and time, and aid in interpreting environmental controls on species distributions [51][52][53][54][55][56][57][58][59][60]. HCA is a classification tool that iteratively partitions samples into groups based on differences in their fossil assemblages. Initially, the HCA algorithm links the two samples with the shortest distances into a group. Next, a new distance matrix is computed from the remaining samples. The sample with the shortest distance to the first group is linked to it. This procedure repeats until all samples are combined into groups and the groups are fused into clusters [50,61]. In this way, each cluster represents a set of samples with similar palynological compositions. The most similar samples will have the lowest Euclidean distances. DCA [62] is a popular ordination technique for detecting gradients of ecological change and relating this variability to underlying environmental factors [63][64][65][66][67]. The relative position of samples within DCA space reflects their biotic similarity. Samples that plot close to one another have more similar biotic compositions than samples that plot far from one another. DCA captures the primary biotic variation along DCA axis 1. Subsequent axes explain smaller amounts of variation. This quantitative approach to biofacies definition permits the fossil data to reveal significant relationships that tell their own story unconstrained by implications from other data. The resulting biofacies are direct products of the statistical analyses, based upon the degree of dissimilarity in fossil composition among samples as measured by the Euclidean distance coefficient. An advantage of this strategy is that the interpretation of external controls on biotic variability is relatively straightforward and achieved through overlaying environmental information onto the cluster dendrogram and ordination plot [47]. A link between biotic patterns and environmental controls is established when the environmental data maps convincingly onto the biofacies interpretations. If there is not a good match between the interpreted biofacies and environmental data, then, the environmental data likely had little influence over biofacies composition. We coded the samples in the ordination by locality, cluster membership, time horizon, paleosol type, and depositional environment to aid in interpreting controls on biotic variability. A second advantage of this approach is that samples and taxa can be plotted together within the same ordination space. Samples that plot close to a particular taxon typically have the greatest abundances of that taxon [47]. This makes it easy to visualize the taxa that characterize each biofacies, and to interpret gradients in biotic composition that can ultimately be related to environmental gradients. All multivariate analyses were performed using the R environment for statistical computing [68]. HCA was performed using the AGNES function from the CLUSTER package [69]. DCA was performed using the DECORANA function from the VEGAN package.
Analytic rarefaction [70][71][72][73][74][75] was used to compare taxonomic diversity (e.g., richness) among the biofacies, localities, paleosol horizons, and depositional environments studied. Rarefaction computes estimates of taxonomic richness and 95% confidence intervals at a standardized, scaled down sampling effort so that comparisons can be made among samples of different sizes. Rarefaction was performed using the program Analytic Rarefaction version 1.3 [76]. In this study, sampling effort is defined by the number of fossil individuals contained within each pooled sample grouping for comparisons among biofacies, localities, paleosol horizons, or depositional environments.

Hierarchical Agglomerative Cluster Analysis (HCA)
• Five clusters, referred to as biofacies A-E are interpreted in the cluster dendrogram (see Figure 4). • A significant branch point at a Euclidean distance of 0.25 separates biofacies A and B from biofacies C, D, and E ( Figure 4). This branch reflects a major break in biotic composition, from the fern and moss dominated samples of biofacies A and B to the brackish and freshwater algae dominated assemblages of biofacies C, D, and E.

•
In general, clusters tend to differentiate samples among the localities and the depositional environments from which they were collected, although overlap exists. The clusters do not cleanly segregate samples of different paleosol types or from different paleosol horizons, although loose groupings are observed (see Figure 4).

•
Biofacies A mainly comprises swamp and lake margin samples from the P3 through P6 paleosol horizons of the Sentinel Hill and Kikiakrorak River Mouth localities. Fern and moss spores dominate, especially Psilatriletes, and comprise 56% of the biofacies. Brackish and freshwater algae, including Sigmapollis, are common and comprise 19% of the total counts in the biofacies (see Figure 4 and Table 2). • Biofacies B mainly contains samples from overbank facies of the P2, P3, P7, and P9 paleosol horizons from all four localities. Similar to biofacies A, biofacies B is also dominated by fern and moss spores (43%) and algae (25%). Unlike biofacies A, Psilatriletes is rarely encountered. Instead, the spore Laevigatosporites is the most abundant genus (15%) (see Figure 4 and Table 2). • Biofacies C mainly contains samples from swamp margins from the P6 and P9 paleosols of the Liscomb Bonebed and Kikiakrorak River Mouth localities. Brackish and freshwater algae (39%) and exotic pollen (26%) are dominant. Diagnostic taxa include Sigmapollis and the pollen genus Aquilapollenites (see Figure 4 and Table 2). • Biofacies D is characterized by samples from the undifferentiated lower delta plain from the P8 and P9 paleosols of the Liscomb Bonebed. Biofacies D contains a high abundance of brackish and freshwater algae (50%) and fern and moss (22%) genera.
Sigmapollis is dominant, composing nearly 40% of the biofacies (see Figure 4 and Table 2). • Biofacies E comprises two samples from overbank and swamp margin deposits of the P1 and P2 paleosols from the North Kikak-Tegoseak locality. It contains the highest proportion of algae genera (70%) observed. Unlike the other algae dominated samples from biofacies C and D, in biofacies E Botryococcus algae, not Sigmapollis, is diagnostic (see Figure 4 and Table 2). samples from biofacies C and D, in biofacies E Botryococcus algae, not Sigmapollis, is diagnostic (see Figure 4 and Table 2).

Detrended Correspondence Analysis
• Coding the samples by biofacies membership reveals that the DCA results largely support the results of the cluster analysis. • Samples from each individual biofacies (A through E) tend to plot cohesively in ordination space, although some overlap exists (see Figure 5A). This suggests that biofacies share aspects of their biotic compositions, as show in Table 2.

•
The major segregation between fern and moss dominated biofacies (A and B) and algae dominated biofacies (C, D, and E) is clearly observable (see Figure 5A). Algae dominated samples have low axis 2 scores, while fern and moss dominated samples have intermediate to high axis 2 scores. The DCA scores of select algae genera (Sigmapollis and Botryococcus) and fern and moss genera (Psilatriletes, Aquilapollenites, Deltoidospora, and Laevigatosporites), as well as the average scores of all taxa from these two respective ecological groups, support the separation of algae from fern and moss dominated assemblages. • One algae dominated sample from biofacies E plots outside the main cloud of algae dominated samples. Its high axis 1 score is driven by elevated abundances of the algae Botryococcus, a genus that is less common in other algae dominated samples (see Figure 5A and Table 2).
When coded by locality, samples from the Liscomb Bonebed, Sentinel Hill, and Kikiakrorak River Mouth plot with low axis 1 scores. These samples largely separate from the North Kikak-Tegoseak samples, which have intermediate to high axis 1 scores. Along axis 2, samples from the Liscomb Bonebed and North Kikak-Tegoseak localities display low scores and separate from the Sentinel Hill and Kikiakrorak River Mouth localities, which plot with higher scores (see Figure 5B). Variations in the abundances of ecological groups and individual taxa among localities are shown in Table 4. Geosciences 2021, 11, x FOR PEER REVIEW 10 of 22  • Coding the samples by paleosol horizon reveals no easily generalized temporal trend, although samples from each individual horizon tend to plot closely in space (see Figure 5C). • When samples are coded by paleosol taxonomy, a weak separation is observed between samples from acid sulfate soils (horizons P8 and P9 from the Liscomb Bonebed) and the aquic entisols and inceptisols that characterize all other samples (see Figure 5D).  Figure 5E). • Swamp margin, lake margin, and overbank facies are dominated by fern and moss genera, which comprise nearly 40% of each environment's biota (see Table 4). Swamp and lake margin samples share many common and abundant genera, including Psilatriletes and Sigmapollis. Variation in the abundances of the lowland tree/shrub Taxodiaceaepollenites and the algae Pediastrum differentiate the swamp and lake margin (see Table 3).
Plotting the scores of taxa within DCA space corroborates these compositional trends.
The average score of fern and moss taxa plots closely to swamp and lake margin samples and nearby to overbank samples, indicating they are common elements of these environments. • Overbank samples are dominated by Sigmapollis and Laevigatosporites; Psilatriletes is rare. The overbank contains a higher proportion of hinterland conifer pollen, reworked dinocysts (~1%), and marine acritarchs and peridinoids (~1%) than the lake and swamp margins (see Table 3).

•
The undifferentiated lower delta plain samples are dominated primarily by brackish and freshwater algae and secondarily by fern and moss genera. Sigmapollis is the most abundant taxon at 34% abundance (see Table 3). This environment contains the highest proportions of in situ and reworked marine elements (dinocysts and acritarchs) observed in the study (~3%), and these tend to plot close to undifferentiated lower delta plain samples in ordination space (see Figure 5E).

Analytic Rarefaction Analysis
• Rarefaction curves displaying richness at different sampling efforts are shown in Figure 6A-D. • Biofacies B has a significantly higher taxonomic richness (61.4 taxa) as compared with other biofacies at the lowest common sampling effort (420 individuals). This is largely due to increased diversity of fern and moss and lowland tree and shrub taxa (see Table 5). Biofacies A, C, and D have lower and overlapping richness values (47.7, 45.1, and 49.1 taxa, respectively). Biofacies E has the lowest richness value (36.5 taxa) overall (see Figure 6A).  Figure 6B and Table 6).  Figure 6C and Table 7). abundant taxon at 34% abundance (see Table 4). This environment contains the highest proportions of in situ and reworked marine elements (dinocysts and acritarchs) observed in the study (~3%), and these tend to plot close to undifferentiated lower delta plain samples in ordination space (see Figure 5E).

Analytic Rarefaction Analysis
• Rarefaction curves displaying richness at different sampling efforts are shown in Figure 6A-D.
• Biofacies B has a significantly higher taxonomic richness (61.4 taxa) as compared with other biofacies at the lowest common sampling effort (420 individuals). This is largely due to increased diversity of fern and moss and lowland tree and shrub taxa (see Table 5). Biofacies A, C, and D have lower and overlapping richness values (47.7, 45.1, and 49.1 taxa, respectively). Biofacies E has the lowest richness value (36.5 taxa) overall (see Figure 6A).

•
The Liscomb Bonebed, North Kikak-Tegoseak, and Sentinel Hill localities display statistically indistinguishable richness values (66.3, 62.9, and 63.2 taxa, respectively) at the lowest common sampling effort (720 individuals). The Kikiakrorak River Mouth locality has significantly lower richness (57.2 taxa) (see Figure 6B and Table  6).  Figure 6C and Table 7).   • Overbank environments have the richest assemblage (57.8 taxa) of any environment at the lowest common sampling effort (380 individuals). This is attributed to greater diversity of algae, fern and moss, hinterland conifer, lowland tree and shrub, and reworked marine taxa. The undifferentiated lower delta plain and swamp margin environments display overlapping and intermediate levels of richness (57.7 and 55.3 taxa, respectively); the lake margin has the lowest sampled richness (37.9 taxa) (see Figure 6D and Table 8).

Environmental Controls on Biofacies Composition and Paleosol Development through Space
Prince Creek biofacies vary along gradients that reflect moisture level (drainage) and marine influence (see Figure 5F). Biofacies A, C, and D occur within the predominantly saturated soils of lake margin, swamp margin, and undifferentiated lower delta plain environments, while biofacies B and E tend to occur in overbank soils that, while still wet, were subject to periodically drier conditions [17,22]. In general, samples from the intermittently drier extreme of the gradient record significantly higher abundances of bisaccate pollen from hinterland conifers. Bisaccate pollen are transported from the uplands and delivered by the river system to the raised and intermittently drier, proximal levees of the overbank setting. Independent evidence for drier overbank conditions come from common occurrences of ferruginous features and Mn oxides, higher bioturbation intensities (reflective of stable soils), the absence of aggregated, zoned soil structures (peds), and increasing Fe/Al geochemical ratios indicative of oxidized and better drained soils [22]. Alternatively, samples from the predominantly wet end of the gradient, which include swamp and lake margins of the lower delta plain, have markedly reduced abundances of hinterland conifer pollen. The lower delta plain is situated in a topographically lower position than overbank environments more proximal to channels, and receives proportionately less bisaccate pollen from the uplands. As noted above, lake margin, swamp margin, and undifferentiated lower delta plain environments tend to be dominated by fern and moss (Psilatriletes), and algae (Sigmapollis and Botryococcus) assemblages and are indicative of wetter habitat conditions that are possibly dystrophic where Botryococcus dominates [77]. Marine taxa, including in situ peridinoid dinocysts and acritarchs, attain their highest abundances in the undifferentiated lower delta plain, lending additional support for the interpretation of wetter conditions along this portion of the gradient and proximity to the coastline. Furthermore, the presence of horizons with high total organic carbon (TOC), pyrite, gypsum, jarosite, Fe-depletion features, drab-soil colors, carbonaceous plant fragments, and a lack of bioturbation in these soils are suggestive of more poorly drained conditions near the coast [17,22].
The gradient along DCA axis 1 also roughly reflects an associated change from the distal (lowland) coastal plain, represented by the low axis 1 scores of Liscomb Bonebed samples, to the more proximal (upland) position of samples from North Kikak-Tegoseak. Indeed, evidence from previous stratigraphic and paleontologic studies suggest that the North Kikak-Tegoseak locality is situated in the most updip position in the study, in a more proximal position relative to the Brooks Range orogenic belt and the fluvial systems that delivered sediment to the distal delta plain, while the Liscomb Bonebed is located in one of the most distal locations along the lower coastal plain or delta plain near the non-marine to shallow-marine transition zone [13,15,[17][18][19]21,22]. Lowland and more distal coastal plain or delta plain localities would be subject to prolonged waterlogging, possibly due to lower topographic relief, a higher water table, seasonal river flooding, annual changes in water table position, and marine transgressions [13,17,21,22,78,79]. Topographically higher positions in more proximal parts of the coastal plain likely experienced better drainage and remained drier for relatively longer periods of time [13,22,78,80]. This interpretation of paleosol-type relationship to topographic gradient corroborates findings from the analysis of stable oxygen isotopes in dinosaur tooth enamel [14]. Suarez et al. [14] concluded that enamel from Pachyrhinosaurus fossils of the North Kikak-Tegoseak dinosaur bonebed locality were enriched in δ 18 O because they foraged on enriched upland conifers. Alternatively, enamel from Edmontasaurus dinosaur fossils of the Liscomb Bonebed was depleted in δ 18 O because Edmontasaurs consumed isotopically depleted plants from along the distal coastal plain.
Biofacies change along axis 2 is controlled by marine influence. The increased presence of in situ marine and brackish palynomorphs, ostracodes [34], Nucula clams [18], as well as pyrite, gypsum, and jarosite [22] suggest that samples with lower axis 1 scores were at times influenced by the input of brackish/marine groundwater near the shoreline. Pyrite formation occurs in soils due to the interaction of iron with sulfate within pore waters [81]. Sulfate commonly occurs within marine, rather than freshwater settings [82]. Jarosite and gypsum are both oxidation products of pyrite [82,83], again suggesting increasing marine influence in paleosols that contain these minerals in abundance. The relative absence of marine and brackish taxa from samples with higher axis 2 scores, coupled with the presence of sphaerosiderite, suggests these samples experienced dominantly fresh ground waters. Unlike pyrite, sphaerosiderite tends to precipitate only when dissolved sulfur concentrations are low and pore waters are fresh [84][85][86][87]. Stratigraphic studies of the Prince Creek Formation at Sling Point, the Liscomb Bonebed, Ocean Point, and nearby localities suggest that this stratigraphy records increasing marine conditions up section, prior to apparent transgression recorded in shallow-marine deposits of the overlying tongue of the Schrader Bluff Formation evidenced at Ocean Point [15,19,21,23].

Environmental Controls on Biofacies Composition and Paleosol Development through Time
Within the confines of the current sampling effort, it is challenging to characterize and quantitatively compare biotic gradients among the studied paleosol horizons. As shown in Table 1, the sampling of depositional environments within each horizon is largely incomplete. In fact, there is no single horizon from which all four depositional environments have been sampled. Likewise, each horizon is represented by only a small number of samples that may not be adequate to estimate or compare the constituent biotas of each paleosol individually. Finally, although the relative ages and updip-downdip relationships of the studied paleosols are relatively well understood [13,15,[17][18][19]21,[27][28][29][30][31][32]34,88], the correlation of paleosol horizons over broad areas is extremely difficult since the extent of any one paleosol horizon is limited due to truncation by channel deposits, the semicontinuous nature of outcrop exposures, and lack of well and seismic control [22]. Despite these challenges to accurate, statistically sound temporal comparisons, the ordination of samples coded by paleosol horizon (Figure 5C) corroborates the interpretation of a transition toward greater marine influence in younger Prince Creek stratigraphy exposed near the site of the Liscomb Bonebed [15,21,23]. Samples from the youngest sampled horizons, P8 and P9, tend to separate from the older P2 through P7 horizons.

Future Research Directions
The multivariate statistical approach advocated in this study can be used to address many remaining questions about the nature of the Prince Creek ecosystem, how it varies geographically and temporally, and whether this variation is typical of coastal plain settings in general:

1.
How much biotic turnover is typical among depositional environments within a paleosol horizon? Does the level of between habitat variability change through time? 2.
Does each paleosol horizon contain a unique palynological/microbiotic signature, or does the biota tend to recur with only slight variations through time? 3.
What does any observed biotic variability tell us about the evolution of the coastal plain ecosystem and if/how physical environmental factors vary though time? 4.
How does the degree of biotic and environmental change observed compare to that in lower latitude Maastrichtian settings, where high frequency, seasonal changes may not be as evident? 5.
How similar are patterns of palynofacies variability and interpreted environmental controls in other marginal marine settings outside of the PCF?
To address these questions, future field research could focus on collecting closely spaced replicate samples from each PCF soil horizon to permit nested statistical comparisons of soils at multiple scales. In this way variability can be examined among replicate samples within each horizon, environment, or locality. A detailed understanding of withinhorizon variability is crucial so that a statistical baseline of change can be established, and then used to evaluate the significance of variability among horizons, environments, or localities [58,65,[89][90][91][92]. With greater sampling intensity, additional multivariate techniques including analysis of similarity (ANOSIM) [56] can be used to test statistically for differences within and among multiple test levels to address the above questions. In addition to field sampling, efforts can be made to compile previously published quantitative palynomorph data from the PCF, for example, [93], to expand the scope of future studies. Finally, we suggest that more quantitative studies of palynofacies in coastal plain ecosystems are needed to better understand whether the variability we observed is typical of these marginal marine settings. The answers to the above questions can be integrated with existing observations from stratigraphy, sedimentology, paleopedology, and geochemistry to provide a more highly resolved view of the Prince Creek ecosystem in Alaska, marginal marine systems elsewhere, and establish well-supported links between environmental and biotic variability.

Additional Uses of Quantiative Biofacies Analysis/Multivariate Statistical Tools
This quantitative approach to biofacies analyses can be used for other purposes, as well as in stratigraphic intervals outside of the PCF of Alaska. Because stratigraphic architecture and environmental change affect fossil assemblages in predictable ways [37,40,47], a biofacies analysis with HCA, DCA, or other ordination techniques provides a useful tool for building interpretations of stratigraphic and environmental architecture [46,48,60] and for regional and intraregional correlation of horizons [64] that are independent of lithological, geochemical, or other data. Although a quantitative biofacies analysis tends to be more common in academic studies, it can also prove useful in building predictive stratigraphic, depositional, and reservoir models for industry purposes [94].
Multivariate statistical analyses can be applied broadly whenever one seeks to summarize quantitative multivariate data, classify groups based on shared similarities of properties, or relate and display statistical relationships among multiple objects. Due to the advent of "big data", tools such as cluster analysis, ordination, and others are increasingly used by geologists to extract patterns from subsurface data. Multiple examples are published that provide illustrative cases. For example, in areas where regional correlation is challenging due to a lack of biostratigraphic data, surface exposures, or seismic data, cluster and ordination analyses can be used to develop chemostratigraphic correlations based on similarities in geochemical, elemental, and isotopic signatures [95,96]. These tools are also useful for analyzing biomarker and other geochemical data to characterize oil families and understand regional differences in petroleum systems [97,98]. Geophysicists are turning to principal component analysis (PCA) and artificial neural networks to evaluate which combinations of attributes extracted from 3D seismic data best reflect hydrocarbon bearing reservoirs [99]. Additionally, development geologists and engineers use multivariate and artificial intelligence tools to understand which reservoir properties are most important in driving both production performance [100,101] and variability across hydrocarbon producing trends.

Conclusions
Cluster and ordination analyses reveal that palynomorph and microbiota of the PCF coastal plain can be categorized into two main assemblage types: (1) fern and moss dominated biofacies characterized by the typically water-logged lake margin, swamp margin, and lower delta plain paleosols, and (2) algae dominated biofacies comprising periodically drier overbank paleosols. Biofacies are arrayed along environmental gradients reflecting moisture level (degree/frequency of water-logged conditions) and marine influence. These findings broadly support previously published depositional models of the Prince Creek coastal plain based on paleosol pedology, geochemistry, and stratigraphy [22]. One of the strengths of this study, and a key difference from previous work, is the application of a robust statistical methodology to define biofacies, visualize biotic variation within and between depositional environments, and link biotic change to habitat change. Comparisons of within-habitat biotic variability suggest that the Prince Creek coastal plain was a dynamic ecosystem subject to frequent fluctuations in environmental conditions. These fluctuations are not only manifested in the sedimentologic and geochemical record, but also reflected in the variable biofacies compositions recorded within and between coastal plain depositional environments. A quantitative biofacies analysis provides an additional independent tool for interpreting changes in ecosystems through space and time, and understanding the stratigraphic, environmental, and evolutionary processes effecting these changes.