The fossil communities used in this study (Figure 1) come from the Pliocene Beaufort Formation and Beaufort age-equivalent “High Terrace Sediments” of the CAA, which comprise unlithified sand and gravel containing organic deposits from braided rivers (Fyles, 1990). These rivers were likely surrounded by forests that may have extended all the way to the Arctic Ocean during parts of the Pliocene, nearly eliminating the tundra biome (Salzmann et al., 2008). Found within the organic deposits are remarkably preserved woody debris from the forests (Devaney, 1991; Fyles et al., 1994; Murphy, 2006; Murphy et al., 2007) up to half a meter thick (Davies et al., 2014), and layers of well-preserved leaf litter comprising leaves, twigs, mosses, pollen, charcoal, and invertebrate remains (Fyles, 1990; Matthews and Ovenden, 1990; Elias and Matthews, 2002; Tedford and Harington, 2003; Rybczynski et al., 2013; Fletcher et al., 2016). The deposits also preserve a fascinating mammalian fauna (Tedford and Harington, 2003; Rybczynski et al., 2013), some of which challenge previous hypotheses for the evolutionary and ecological history of modern lineages (Rybczynski et al., 2013). Thus, the study of the climate and environment of these fossil communities may also improve our understanding of the past ecology that shaped the evolution of modern mammalian forms.

Dating of the Beaufort Formation and “High Terrace Sediments” is complicated by the combination of unconsolidated, often coarse sands, proximity to the pole, and absence of volcanic sediments in the formations. These sites have been dated using techniques including biostratigraphy of the flora (Matthews and Ovenden, 1990), and fauna (McNeil, 1990; Fyles et al., 1994; Tedford and Harington, 2003) the relative stratigraphic position (Fyles et al., 1994) and timing of eustatic highs (Brigham-Grette and Carter, 1992), magnetostratigraphy (Fyles et al., 1991), Strontium (Sr) isotope analysis (Kaufman et al., 1990), amino acid racemization (Brigham-Grette et al., 1987; Fyles et al., 1991), and more recently absolute dates obtained from terrestrial cosmogenic nuclides (TCN; Rybczynski et al., 2013). Justification for our current interpretation of the ages of the sites in this study follows below, and is shown in Figure 2, however, the errors on these estimates are substantial and a more tightly constrained chronology of the sites would increase their value to Pliocene climate modelers and for understanding the geologic and evolutionary history of the CAA.

Both Fyles Leaf Beds (FLB) and Beaver Pond (BP) sites from the “High Terrace Sediments” in this study, have been absolutely dated using TCN (Rybczynski et al., 2013), yielding minimum ages of 3.8 +1/−0.7 and 3.4 +0.6/−0.4, respectively, and thus overlap within one standard error. The relative elevation between the significant peat layers at the two sites, both at ~400 m elevation, supports an approximate coeval deposition. In addition, they share specific informative floral elements such as Paliuris and shared components of faunal assemblages (Rybczynski et al., 2013).

The date for the Beaufort Formation at Meighen Island was previously considered to be ~3 Ma based on the estimated age of the marine incursion layer that separates the lower sediments of the Beaufort Formation on Meighen Island (MI-1) from the upper sediments of Meighen Island (MI-2), and this age is still often cited in the literature (e.g., Elias and Matthews, 2002; Haywood et al., 2016). This estimated age was primarily based on the discovery of a mixed molluscan fauna, comprising both Atlantic and Pacific Ocean species, suggesting a date of deposition proximal to the opening of the Bering Strait, secondarily evaluated using magnetostratigraphy (Fyles et al., 1991). It was compared to the then best known site with a similar mixed flora and magnetostratigraphic record—the Tjornes Beds, Iceland. The dates of the Tjornes Beds have been revised several times since (Fyles et al., 1991). The most recent, and perhaps complete dating of the section, places the mixed molluscan flora recognized at this site, at the base of the Serripes Zone at 4.4–4.5 Ma based on dinoflagellate cyst stratigraphy and paleomagnetic correlation with the top of the Nunivak subchron (Verhoeven et al., 2011).

The opening of the Bering Strait and ages of mixing of the Atlantic and Pacific molluscan faunas, has also been re-evaluated (Marincovich, 2000; Marincovich and Gladenkov, 2001; Gladenkov et al., 2002). Marincovich (2000) suggested that there were multiple mixing events, with one following the first opening of the Bering Strait (4.8–5.5 Ma), and a second mixing event due to the change of direction of flow through the strait, arguably due to shoaling of the Central American Seaway (~3.6 Ma; Although see Ramírez et al., 2016). These two periods also coincide with eustatic episodes (Marincovich, 2000). The thickness of the marine beds (<100 m) was considered to suggest a short period of deposition (Fyles et al., 1991), and this may be a “better fit” for the shorter high stand at ~3.6 Ma. However, the sequence of opening, shoaling and complete closure of the key seaways is not agreed on in the literature, and the interchange of fauna may have persisted beyond initial mixing.

The ages of the Meighen Island marine incursion layer have also been subject to other dating attempts. Sr isotope dating of marine shell estimated an age range of 2.5–5.1 Ma (Kaufman et al., 1990), and, amino acid racemization indicated the ground temperature following deposition must have been freezing or below to promote the degree of amino acid preservation observed thus suggesting the shells were not subject to high temperatures, such as those of the Early Pliocene, for an extended time (Brigham-Grette et al., 1987; Brigham-Grette and Carter, 1992). The ages were also constrained for the marine incursion by biostratigraphic correlation of marine foraminifera to between 2.4 and 5.33 Ma (McNeil, 1990). Together, these lines of evidence are consistent with deposition of the marine incursion layer on Meighen Island at approximately the boundary between the Zanclean and Piacenzian. Thus we consider the MI-1 flora to be ~3.6 Ma or older, and MI-2 ~3.6 or younger.

The early reports of the Banks Island Beaufort Formation (BI) age was suggested relative to Meighen Island, while noting that some of the flora indicated an Early Pliocene age, and that the composition more generally suggested an age younger than Miocene (Fyles et al., 1994). Having reviewed the dating of the MI marine incursion we consider that this may also impact the likely age of the Beaufort Formation at BI. Thus, the BI Beaufort Formation age remains poorly known.

The oldest and youngest ages for all sites in the study are additionally constrained by key elements of the flora. For example, five-needle Pinus have been recorded in the Banks Island flora, both Meighen Island floras and the BP flora and is considered an indicator of pre-Quaternary environments in the High Arctic. The presence of the myriaceous taxa Myrica (Gale) and Comptonia also strongly suggests a Pliocene or older age, although their precise Last Appearance Datum (LAD) is not known. Epiprenum crassum's LAD in Europe is Early Pliocene (Matthews et al., 2003), thus the presence of this taxon at BI and MI-1 but not MI-2 despite significant sampling and investigation, while acknowledging that absence of evidence is not evidence of absence, may lend support to the Meighen Island marine incursion occurring at the Zanclean-Piacenzian boundary.

Faunal elements lend another line of evidence constraining the ages presented here. Tedford and Harington (2003) compared the vertebrate elements of BP to the Early Pliocene Yushe Basin faunas of northeastern China, and suggested an age of 4–5 Ma, as well as an active interchange across the Bering Strait. The insect fauna generally does not have biostratigraphic significance, but where some constraints are suggested, they confirm the Pliocene ages assigned to the localities (Elias and Matthews, 2002).

Based on these considerations, we suggest that the ages of these sites may fall at approximately the Zanclean-Piacenzian boundary, however, due to substantial uncertainties, only the two Meighen Island sites can be assigned a definite sequence relative to each other.

We adopted the plant taxon lists for this study (Data Sheet S1), from Fyles et al. (1994), Matthews and Ovenden (1990), Matthews and Fyles (2000), and combined these with previously unpublished taxa lists. These identifications were derived from a wide range of remains of plant macrofossils preserved in sediments best described as peats through to leaf beds with little decomposition evident, protected by permafrost that has covered the area for much of the period since their deposition. Each locality herein may refer to one site or a group of sites in the same region that is suspected on stratigraphic evidence to be of a similar age, and is considered by Matthews to be one flora.

The identification of many of the elements of the floras to extant taxa provides the opportunity to apply both new methods and established methods to estimate the paleoclimate. The recently devised Climate Reconstruction Analysis using Coexistence Likelihood Estimation (CRACLE) derives the most likely climate for the combination of floral elements identified using joint likelihood functions, based on modern distribution data from the GBIF database (Harbert and Nixon, 2015). The CRACLE estimate median differs from the modern WorldClim dataset by an average of <1.5°C, significantly improving on the error from previous paleovegetation methods (Harbert and Nixon, 2015), while integrating more of the available community data thus mitigating the effect of mixed floras, and addressing some methodological concerns of the commonly implemented coexistence approach (CA; For detailed discussion see Grimm and Potts, 2016). The CRACLE estimates are compared to the CA (Full results see Data Sheet S2; Mosbrugger and Utescher, 1997), which calculates climate estimates based on the total overlap of climate tolerances of the nearest living relatives of the flora non-probabilistically, but modified from the classical approach in that we did not use nearest living analogs for extinct species, we only used those that are extant (see Data Sheet S1).

For our climate reconstructions, three subsets of the data were derived from the plant taxa listed in Data Sheet S1, which capture a trade-off between increasing sample size and decreasing precision in identification (Table 1). The three subsets of data include: (1) confident identifications of extant taxa at the species level (referred to from this point as species), (2) confident identifications of extant taxa at the species level and taxa that compare well to extant taxa at the species level (referred to from this point as species cf.), and (3) those that are confidently assigned to the genus level with extant members (genera). The occurrence records for each list were accessed from the Global Biodiversity Information Facility (GBIF; Hijmans et al., 2005) via Darwin Core Archive. Those observations for which the value for “hasCoordinate” was false and/or “hasGeospatialIssues” was true, were removed. The remaining occurrence records were used for the input into the CRACLE code (Harbert and Nixon, 2015) for R (R Core Team, 2013), with some modifications to the body of the script guided by Harbert (R library in development).

The output for the provided CRACLE R-code includes the results of the parametric and non-parametric joint likelihoods (P-CRACLE and N-CRACLE), the ranges derived by “strict” CA (100% only), and the ranges derived from a “relaxed” (90%) CA (Mosbrugger and Utescher, 1997) using the GBIF database observations and WorldClim as inputs for the 19 climate variables used by BIOCLIM (http://www.worldclim.org/bioclim; see Data Sheet S2 for full results). Of these climate variables, six climate variables from Harbert and Nixon (2015) verified against WorldClim were assessed in detail: Mean Annual Temperature (MAT), maximum temperature of warmest month (MAX), minimum temperature of the coldest month (MIN), Mean Annual Precipitation (MAP), precipitation of the wettest quarter (WET), and precipitation of the driest quarter (DRY). As per the original tests of the CA against fossil flora (Mosbrugger and Utescher, 1997), using all of the coexisting taxa produced non-overlapping climate ranges. For this reason the strict (all data from all taxa) CA was not used, and further discussion refers to the CA method as proposed by Mosbrugger and Utescher (1997), which calculates the envelope in which the most taxa are able to coexist (referred to as the relaxed CA in code output). In addition, we investigated if removing non-woody taxa, as per Harbert and Nixon (2015), differed from our complete taxa list results by running CRACLE by woody species only.

In order to investigate whether the dissimilarity in community structure depends upon a dissimilarity in climate, we tested the relationship between climatic results and dissimilarity indices for the five sites. The lists of plant taxa identified across the CAA facilitate an investigation of the shared (beta) diversity between sites, and thus the similarity or dissimilarity of the floras. The use of dissimilarity indices suitable for paleontological applications, in this instance the Raup-Crick Dissimilarity Index, a weighted measure of community dissimilarity using Monte Carlo simulations (Raup and Crick, 1979; Chase et al., 2011), provide evidence with regard to the interconnectedness of the communities studied and the conditions in which they formed.

The Raup-Crick analysis was conducted separately on vascular plants at the level of species and species cf. for correlation with climate variables as was used for climate analysis (1:0, sites as rows, Data Sheet S3). The implementation follows Chase et al. (2011). Dendrograms of the Raup-Crick indices were generated using h-clust with agglomerative clustering using complete linkage (furthest-neighbor joining), and the Unweighted Pair Group Method with Arithmetic Averaging (UPGMA) in vegan (2.4-1) for R (Oksanen et al., 2015). Complete linkage is a conservative method, which tends to underestimate similarity between recognized clusters, whilst UPGMA is the more commonly applied method in paleoecology. The resulting patterns of clustering were then compared.

The difference in the six climate variables, pair-wise by site, was tested against the dissimilarity index results using the non-parametric Kendall's Tau. This approach was selected as a Shapiro–Wilk test resulted in a large number of non-normal distributions of variables and it produces more accurate P-values for small samples. Kendall's Tau was calculated using “corr.test” in psych (1.6.9) for R (Revelle, 2016). The species derived climate variables, and formed pairwise combinations, were tested against the species dissimilarity indices. The species cf. derived climate variables, and formed pairwise combinations, were tested against the species cf. derived dissimilarity indices. To understand if the results could be due to low statistical power, the climate estimates from species were examined for clustering by recursive partitioning using “rpart” (4.1.10) for R (Therneau et al., 2015) and compared to the clustering of dissimilarity indices from species. The basis of node branching was compared with results from a principal components analysis.

The warmest of all sites investigated in this study was BI, with the contribution to these higher temperatures relative to the other floras, more from warmer winter temperatures, than from warmer summer temperatures i.e., temperature of the coldest month was much warmer at BI, whereas temperature of the warmest month was only slightly warmer than the other sites (Figure 3, Data Sheet S2). This pattern of reduced thermal amplitude, with increased MAT is consistent with previous studies, and is evidence for the internal consistency of this method (e.g., Hernández Fernández, 2006; Hernández Fernández et al., 2007). The next warmest were MI-1, FLB, and MI-2, whilst BP had the coolest MAT estimates. WET and DRY were fairly consistent between sites and WET tended to approximate or be slightly higher than three times the precipitation of DRY. The variability in the MAP results between species and species cf. for BI makes it difficult to compare the sites for this variable. The range of estimates suggest that this site is either approximately the same as MI-2, BP and FLB or considerably wetter.

Between the analyses conducted using species only, or species cf. identifications, DRY, WET, MAX, and MIN were consistent. For MAP MI-1, MI-2, and BP were similar, but FLB was drier when using species cf. and precipitation at BP was nearly twice as high if species cf. is used. MAT was warmer with species cf. for BI ~+2°C, MI-2 ~+1°C, and BP ~+0.5°C; MI-1 was very similar between species and species cf. and FLB was slightly cooler at its minimum bound, but barely at the maximum.

For genera level identifications, N-CRACLE results were so broad as to be uninformative (e.g., −30° to +30°C). The genera level results for P-CRACLE and CA were warmer than those derived from the species and species cf. list, which were, in turn warmer than estimates derived from species level identifications. Genera level P-CRACLE tended toward slightly drier results, but with variability between sites and between MAP, WET, and DRY.

The Raup-Crick analyses resulted in different clustering of communities (Figure 4, Table 2). The MI-1 flora clusters with BP and MI-2 with FLB by species, however, the clustering changes when analyzed using species cf. such that MI-1 clusters with FLB and MI-2 with BP. BI is always the most dissimilar locality. Complete linkage clustering and UPGMA yielded similar hierarchical clustering patterns but changed the height of the bracket joining MI-2 and BP to MI-1, and FLB such that complete clustering reported a more positive value indicating a moderate dissimilarity between the two clusters, and UPGMA closer to 0 suggesting the distribution of species between the clusters was not distinguishable from the random draw of the Monte Carlo simulations. Dissimilarity in community assemblage did not correlate by rank with any climate variable or combination of climate variables tested with significance (p < 0.05), even though Banks Island was the most dissimilar in terms of species assemblage and had the highest MAT. Further investigation by examining the clustering of the climate data, showed that the upper and lower bounds of mean annual temperature made the largest contributions to clustering of the sites by climate, and that these clusters were not of the same structure as the clustering of sites by dissimilarity (Data Sheet S4).

The majority of previous climate estimates have focused on a key site for vertebrate fossils, Beaver Pond. The difference between the MAT and MAX estimates for BP and previous estimates (Elias and Matthews, 2002; Ballantyne et al., 2006, 2010; Csank et al., 2011a,b), including those derived using related methods (Ballantyne et al., 2010), are relatively small, often insignificant as they fall within error. Previous MAT estimates from tetraether, 0.6°C ± 5, and paleovegetation, −0.4°C ±4.1 methods are within the range of our estimates by N-CRACLE, for both species 0.8–1.5°C, and species cf. 1.5–1.8°C, and CA, −1.7–5.6°C and 0.5–5.7°C. Ballantyne et al.'s (2010) tree ring isotopes, at −0.5°C ± 1.9, are within the range estimates for species both by N-CRACLE and CA.

Csank et al. (2011a) expanded on Ballantyne et al.'s (2010) study at BP, analyzing ~250 rings and using a mechanistic model to estimate temperature from the isotopic signature. The MAT estimates from Csank et al. (2011a) −1.4°C ± 4.0, include the range of the CRACLE interval, and overlap with the CA. Growing season (June-July) temperature estimates, 15.8°C ± 5.0, were comparable to our species and species cf. MAX interval, 18.2–21.1°C/~20°C, when accounting for the broader window for which the average is calculated. If Csank et al.'s (2011a) and Ballantyne et al.'s (2010) estimates are considered relatively cooler, despite being within error, it may be explained by the uncertainty regarding the isotopic composition of source water leading to potential biases in the these estimates (Sugimoto et al., 2002). For example, a cooler estimate could reflect high use of snow or permafrost meltwater, or magnification of the winter precipitation signal in groundwater. This bias likely has a greater effect on fen environments as posited for BP (Csank et al., 2011a).

The mollusk growth season (May to September) temperature of surface water from Csank et al. (2011b), using conventional and clumped isotopic approaches, resulted in estimates of 14.2°C ± 1.3 and 10.2°C ± 1.4, each with their own caveats. Among these estimates, the upper limits of temperature sampled across the whorl were comparable to our MAX estimates (18.5–21°C). However, isotopic estimates of mollusk growing season temperatures are not necessarily comparable because the paleovegetation records above surface air temperatures, whereas the mollusk shells record surface freshwater temperatures, which are expected to be cooler by at least 1°C (Fricke and Wing, 2004). Air temperature explains only ~38% of maximum summer surface temperature in a modern boreal lake system (Keller, 2007), which potentially explains the relatively cooler temperature estimated by Csank et al. (2011b). In addition to this, and source water composition, other characteristics of the lake, such as water clarity, may also be important for reconciling differences of lake water and surface air temperatures (e.g., Snucins and Gunn, 2000).

Elias and Matthews (2002) used mutual range methods for beetle fauna to estimate MAX and MIN for BI, MI (stratigraphic positions relative to the marine sequence, unknown) and BP. They also found little difference between the sites in terms of MAX (BI 13.5°C, BP 12.4°C, MI 12.9 12.8 12.6°C), however, the estimates are around 7°C lower than presented here. MIN estimates were more variable across the archipelago (BI −21°C, BP −26.6°C, MI −21°C −27.3°C −26.9°C), and were 10–15°C colder than our estimates, but were presented with less confidence than the MAX range (Elias and Matthews, 2002). These estimates point to a clear disagreement between estimates from paleovegetation and isotopic records as compared to paleofaunal mutual range estimates. Warmer estimates from paleovegetation and cooler estimates from the beetle mutual range methods have also been identified at the Late Pliocene Lost Chicken Mine locality in Alaska, at which the presence of Sambucus suggests January temperatures 10–15°C warmer than today, but the beetle mutual range method suggested winter minima only 2.3°C warmer than modern (Matthews et al., 2003). This discrepancy may be due to incomplete knowledge of beetle thermal and distributional ranges given that some species have been observed beyond the official ranges used in the mutual range methods (J. V. Matthews, Jr., Pers. Obs.) or relative insensitivity to minimum temperatures due to behavioral avoidance. This apparently common, if not systematic, difference between floral and faunal estimate warrants further investigation.

The site that is the most distant from the others (BI), was warmest according to N-CRACLE (Figure 3, Data Sheet S2). BI was ~20°C warmer than present (see Table 3 for present values), and the more northerly Meighen Island sites were ~25°C warmer, FLB and BP sites ~23°C and ~20°C warmer than present. The modern difference in MAT between the Eureka Climate Station on Ellesmere Island (−18.8°C), and Sachs Harbour Climate Station on Banks Island (−12.8°C), is 6°C (Environment Canada, 2015). The terrestrial temperature gradient from pole to equator is known to be lower during the Pliocene (Salzmann et al., 2013). Our results suggest the Pliocene temperature gradient may be about half of modern in this region. Under present-day conditions, the difference in temperature between the northerly sites (MI-1, MI-2, BP, and FLB) and BI, is exacerbated by the northern sites' proximity to the thickest Arctic perennial sea ice. Strong cooling of the terrestrial sites occurs directly from the high sea ice, and due to high surface albedo depressing their modern temperatures compared to other locations at the same latitude (e.g., −16.5°C, Vavilov Climate Station, Severnaya Zemlya; Bassford et al., 2006). A lack of or large reduction in permanent sea ice as suggested for the Pliocene (Ballantyne et al., 2013) may have exaggerated the reduction of the temperature gradient between BI and the higher latitude sites. However, when comparing the results from these sites it is necessary to note that although they are all Pliocene, the resolution of dating (see Dating Constraints and Uncertainties) cannot distinguish if they were roughly coeval or represent cooler and warmer periods, perhaps Pliocene glacial and interglacial periods. If so, BP, for instance, may be from a cooler period than its geographically near neighbors.

All sites are currently located in a very dry region (98–130 mm MAP). However, as the geography of the CAA was significantly different during the Pliocene compared to present (Rybczynski et al., 2013), and shorelines advanced and retreated within the time periods represented by these flora (as evidenced by the marine layer that separates MI-1 and MI-2; Fyles et al., 1991), the effect on the precipitation within the region is complex to predict. While, as expected, the analysis estimated precipitation was much higher than today at all sites, the differences in precipitation between the sites are difficult to integrate into our conceptual model of the region given current understanding of the age of the sites and geography of the region at the time of deposition.

The most comprehensive comparison of Pliocene climate models, between each other and with proxy data, is the Pliocene Modelling Intercomparison Project (Hill et al., 2014). Strongly amplified polar warming despite moderate atmospheric CO₂ levels is a notable feature of the time slab studied (Miller et al., 2010; Haywood et al., 2016). Although the models do produce polar amplification, proxy data suggests these are underestimates of the true level of warming (Melles et al., 2012). Indeed, temperature underestimates range from a few to 20°C in northern high latitudes in models such as HadCM3, MIROC, and CCSM4 (Haywood et al., 2010, 2013, 2016). The temperature estimates presented here largely agree with previous estimates (see Comparisons to Previous Estimates), or are warmer than those (Elias and Matthews, 2002). As a result, the addition of this new proxy data will exacerbate the observed proxy-model mismatch for the northern high latitudes. A caveat of comparing our results to the model data outputs is that, although many of our sites' estimates overlap the time slab of these studies (3.264–3.025 Ma), they cannot be precisely dated to this time and thus orbital cycles, variations in atmospheric CO₂, changes in ocean gateways, and many other dynamic processes with large impacts on climate, may generate considerable differences (Feng et al., in press). Additionally, the terrestrial vegetation record for the Pliocene High Arctic may be biased to thermal optima in interglacial periods due to better conditions for higher plant productivity.

It is often necessary in paleoclimatology to use higher levels of identifications where species level is not known, the species is not extant, or may be misidentified as an extant species from limited materials. Given that CRACLE is a recently developed method, and the substitution of genera for extinct species has been criticized for related methods (Grimm and Potts, 2016), we investigated potential concerns in applying the method to different taxa level identifications by extending CRACLE to non-woody species, and using genera level identifications.

Harbert and Nixon (2015) limited their analysis to woody taxa only, but suggest that this is not a necessary limitation on the method. We hypothesized that species level CRACLE estimates would not differ significantly between those based on all vascular plants, and that based on trees and shrubs only. This would be of benefit as not all fossil flora are characterized by abundant and taxonomically rich woody species. Although our hypothesis is largely supported, there are exceptions. The estimates for MAP tend to be somewhat wetter for woody taxa only compared with all species, ranging from a non-discernable difference at FLB to within 50 mm for MI-1, MI-2, and BI, and the largest difference of ~400 mm at BP (Data Sheet S2). Estimates for woody taxa only were also warmer at BP, thus choosing only woody taxa or all vascular species results in a large difference in climate estimates at this site. Overall the difference between woody species only and all species was smaller than the difference between species and species cf. estimates for temperature, but not for precipitation.

The further back in deep time the fossil locality investigated was deposited, the more often identifications can only be made to a genus level, and the ranges of modern species from the same genus have often been used as nearest living analogs. It was predicted that the genus level estimates would be warmer than species level predictions due to the latitudinal biodiversity gradient, biasing diversity within a genera toward lower latitudes. This was supported for P-CRACLE and CA for MAT, MAX, and MIN. The genera level estimations for N-CRACLE are massively wide as to be uninformative e.g., BP MAT. −27°C to 30°C. These results suggests that the use of genera level identifications of fossil taxa are at best biased by global patterns of diversity in the case of the CA and P-CRACLE temperature estimates, if not entirely uninformative as in N-CRACLE, when applied to mutual range methods. Targeted investigation of this effect over a greater latitudinal range would benefit the field of study, but investigators of high latitudes should be especially cautious when interpreting these results.

Comparing precipitation estimates, the N-CRACLE estimates remain uninformatively wide. Within P-CRACLE, DRY estimates are drier for genera level identifications than all other results, and WET are in the same general range. The MAP estimates were mixed. FLB fell within the range of estimates from the other taxa for P-CRACLE, while the other sites were slightly to somewhat drier. This bias appears weaker than the temperature bias, conforming to the results of studies such as Moles et al. (2014) and Butterfield and Munson (2016), which show a stronger association of plant functional traits with temperature than precipitation, if, as posited by phylogenetic niche conservatism (e.g., Donoghue, 2008), these traits are often conserved amongst higher taxonomic rankings.

Harbert and Nixon (2015) suggest CRACLE, and other taxon-based mutual range methods, are most applicable and precise for floras closer to modern (Quaternary) due to sensitivity to extinctions, whilst climates in deep time are more accurately estimated by methods such as leaf foliar physiognomic analyses, which are based on physiological adaptations to climate and do not necessarily require knowledge of phylogenetic relatedness to modern flora (Li et al., 2016). However, foliar physiognomic analyses are complicated by their own biases including taphonomic (e.g., Greenwood, 1992, 2005; Burnham et al., 2001; Kowalski and Dilcher, 2003; Spicer et al., 2011) and methodological (e.g., Jordan, 2011; Milla and Reich, 2011) concerns. The time period from this study lies at the intersection of the projected utility of these two families of methods, thus application of foliar physiognomy to validate CRACLE estimates of climate in the Pliocene may provide valuable insight to disentangle how ecosystem assembly and individual physiognomy respond to climate.

It is evident that the magnitude of the range between the upper and lower boundary for each of the climate estimates given varied between sites and between climate variables within a site (Figure 3). A narrow range estimate of the climatic inferences reflects that there is a narrow range for that climate variable in which those taxa are most likely to coexist. Harbert and Nixon (2015) state that the estimate is influenced most by the most narrowly distributed taxa, and thus the narrower estimates may reflect taxa included at those sites with narrow climate tolerances. However, Harbert and Nixon (2015) did not test the effect of non-overlapping taxa on the range of the climate estimates produced by CRACLE. As all estimates herein have non-overlapping results in the strict implementation of the coexistence approach (See Data Sheet S2), this is an alternative source of influence to be investigated. We found the relationship between the number of taxa at the site or degree of non-overlap of the strict coexistence estimates, and the narrowness of the climate range produced for each variable was unclear. Application of this method to the past would benefit from a dedicated study of the impacts of non-overlapping ranges on climate estimates.

Our analysis suggests that other factors, in addition to climate, may drive the patterns of assemblage of vascular taxa within the CAA during the Early to mid-Pliocene. Within the limits of the range in climate within the region, either the impact of climate is more complex than captured, or other biotic and abiotic factors not measured here contribute a stronger control, or some combination of the two.

Previous studies (e.g., Fyles et al., 1994) have considered the BI flora to represent a boreal forest ecosystem, while growth rings from wood on Meighen Island showed indications of being near tree line. Kuc (1974) and Matthews (1987) also considered the macroflora of Meighen Island indicative of sparse, although diverse (Fyles et al., 1991), conifers at forest-tundra boundary. Matthews and Fyles (2000) describe the BP flora to be largely composed of plants suggestive of an open larch dominated forest-tundra environment containing an extinct species of larch (Larix groenlandii). A more detailed examination of the FLB flora is yet to be published, however it is currently also considered to be located in a forest-tundra ecosystem, featuring the same extinct species of larch (L. groenlandii), abundant leaves of dwarf and shrub birch (Betula spp.) and shrubs of mountain avens (Dryas octopetala) not currently found in the CAA (Matthews and Telka, Pers. Obs.).

Most studies of tree line have focused on regional peculiarities or particular species, rather than generalizable perspectives (Körner and Paulsen, 2004), however at least two global studies have been conducted. In a field study of soil temperature and tree line in modern systems, Körner and Paulsen (2004) found tree line occurs around 6.7°C (±0.8 SD) globally, and that daily means of air and soil temperature are almost the same, between 6 and 7°C. A later study by the same authors (Paulsen and Körner, 2014) developed a model for predicting global tree line, and validated it against satellite observations of tree line position globally. This study also found the mean annual temperature, in this study ~6.4°C, to be a key explanatory variable, along with minimum growth days of 94.

Confidence in applying modern climate tree line boundaries to the past is dependent on an understanding of the causal relationship between temperature and growth form. Many proposed mechanisms for the constant soil temperature—tree line relationship center on soil decomposition rates, low nutrient mineralization rates and their assimilation by trees (See discussion in Scott et al., 1997). This is supported by evidence that photosynthetic rate is not limiting compared with decreased utilization of photosynthetic products in response to low temperatures (Grace et al., 2002). However, these hypotheses alone do not explain the loss of trees specifically, while low shrubs may thrive. Observations have shown that for a given ambient temperature taller flora experience near ambient air temperatures and cooler soil temperatures than do low-statured flora, due to the low flora's aerodynamic decoupling from ambient temperatures above, to form a microclimate around themselves and their soil (Grace, 1988; Scherrer and Körner, 2011). Because trees are not as efficient as low statured plants at generating a warmer microclimate for themselves they are more immediately impacted by lower nutrient acquisition and assimilation rates as a result of low temperature. This suggests a physical and physiological explanation for tree line that should hold if applied to the past.

Combined with the macrofossil interpretations above, this suggests the BI flora was just within this tree line temperature threshold, while our four more northerly sites were outside this range. Thus, the relatively small difference in temperature may be within a critical range to drive the stable and large dissimilarity in community assembly only between BI and the MI-1, MI-2, BP, and FLB cluster, while non-climate influences were stronger drivers of dissimilarity within the MI-1, MI-2, BP, FLB cluster.

The habitat preferences of the species that differ between MI-1, MI-2, BP, and FLB provide evidence for how abiotic differences may have influenced the composition of these communities. For example, taxa at the Meighen Island sites that are known for their role as pioneers include Alnus alnobetula, A. incana, and Betula populifolia, none of which were found at BP, which was likely stable for the duration of peat deposition estimated at 49,000 years (Mitchell et al., 2016). Thus, the stage of the vegetation post-disturbance, and any time averaging of the assemblages may also impact the composition. The dominance of larch (L. groenlandii) at BP and FLB compared to Meighen Island where larch is present but does not dominate, may reflect differences in soils. For example cation, nitrogen and phosphorus availability and water table depth were key factors in a model of evergreen-deciduous-evergreen dominance along mineral supply gradients in boreal peatlands (Givnish, 2002).

Despite these differences, the presence of taxa with habitat preference for wetlands or stream sides is common to all sites. This was confirmed as a Raup-Crick analysis run on only non-wetland species or non-wetland species cf. showed that all sites showed increased dissimilarity except BI, which became either more similar to FLB, or more similar to FLB and MI-1, respectively, although still highly dissimilar according to the index (Data Sheet S3). This presence of taxa with habitat preference for wetlands or stream sides may represent taphonomic bias toward preservation of plant materials in those environments. Occasional exceptions to this pattern, such as the presence of D. octopetala—a taxa more commonly found on dry and rocky sites, suggest transport either of the whole assemblage from multiple sources, deposited in one locality, or transport of the non-wetland components into the wetland system, where they were preserved—an interpretation favored by suitable conditions for preservation within wetland systems.

Thus, we consider dissimilarity to be driven primarily by the temperature threshold for tree line where that is crossed, but by other abiotic factors where the compared sites are below that threshold. Low dissimilarity between the four northern sites is explained by bias toward preservation of wetlands, similarity of climate and the proposed paleogeography of the region at the time featuring a continuous landmass where the archipelago is now incised, thus allowing floral and faunal exchange.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer SF and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.