Assessing the composition of fragmented agglutinated foraminiferal assemblages in ancient sediments: comparison of counting and area-based methods in Famennian samples (Late Devonian)

Benthic foraminifera have been used as proxies for various paleoenvironmental variables such as food availability, carbon flux from surface waters, microhabitats, and indirectly water depth. Estimating assemblage composition based on morphotypes, as opposed to genusor species-level identification, potentially loses important ecological information but opens the way to the study of ancient time periods. However, the ability to accurately constrain benthic foraminiferal assemblages has been questioned when the most abundant foraminifera are fragile agglutinated forms, particularly prone to fragmentation. Here we test an alternate method for accurately estimating the composition of fragmented assemblages. The “cumulated area per morphotype” method is assessed, i.e., the sum of the area of all tests or fragments of a given morphotype in a sample. The percentage of each morphotype is calculated as a portion of the total cumulated area. Percentages of different morphotypes based on counting and cumulated area methods are compared one by one and analyzed using principal component analyses, a co-inertia analysis, and Shannon diversity indices. Morphotype percentages are further compared to an estimate of water depth based on microfacies description. Percentages of the morphotypes are not related to water depth. In all cases, counting and cumulated area methods deliver highly similar results, suggesting that the less time-consuming traditional counting method may provide robust estimates of assemblages. The size of each morphotype may deliver paleobiological information, for instance regarding biomass, but should be considered carefully due to the pervasive issue of fragmentation.


Introduction
Reconstructing paleoenvironments is a key issue for understanding the dynamics of the geosphere and its impact on the biosphere. However, the task becomes increasingly difficult when going back in time. Transfer from modern analogues becomes increasingly problematic with geological time because of the lack of modern equivalents for some extinct groups. Even abiotic proxies such as geochemical data are rarely unambiguous (for instance, oxygen isotope data are affected by temperature of ambient seawater, salinity, and global ice volume) and their interpretation in deep time is subject to debate. The combination of several proxies may be a way to circumvent these problems (Zou et al., 2012;Kumpan et al., 2014;Finkenbinder et al., 2015;Matyja et al., 2015) and enrich our understanding of environmental variations in remote periods of time.
Assemblages of benthic foraminifera have been proposed to trace environmental variations, including nutrient and more generally food availability, oxygenation rate, carbon flux from surface waters, and indirectly water depth (Jones and Charnock, 1985;Loeblich Jr. and Tappan, 1987;Meyn Published by Copernicus Publications on behalf of The Micropalaeontological Society. and Vesperman, 1994;Severin, 1983). In particular, because of the sensitivity of benthic foraminifera to changes in oxygen and organic carbon content in the bottom water, they are prone to respond to relative sea-level variations impacting these parameters (Jorissen et al., 2007). The study of benthic foraminiferal assemblages may thus be pertinent for the characterization of transgressive-regressive patterns in continuous stratigraphic framework (Reolid et al., 2010).
The Late Devonian is such a period, characterized by transgressive-regressive trends and punctuated by major biotic crises (Johnson et al., 1996;Walliser 1996;Haq and Schutter, 2008). Environmental variations throughout the period are complex and are therefore still under debate (Joachimski and Buggisch, 2002;Racki, 2005). Applying an additional proxy, such as benthic foraminiferal assemblages, would thus be of use to better constrain the sequence of environmental changes occurring through the period (Hol-cova, 2004). However, the pertinence of considering benthic foraminiferal assemblages has been questioned. During the Devonian period, the most abundant foraminifera have agglutinated walls (Fig. 1), and the particularly fragile tubular forms dominate in some levels ( Fig. 1a, b). This morphotype is thought to trace deep environments (Nagy et al., 1995;Setoyama et al., 2011). However, because of the fragility of these tubular forms, their estimate in the assemblages may be biased due to fragmentation, therefore artificially biasing their representation in the assemblages (Murray et al., 2011).
We propose here an alternate method to assess the composition of benthic foraminiferal assemblages, circumventing the issue of fragmentation. The composition of the assemblages was quantified based on the cumulative area of each morphotype. We address the following issues. (1) Does the fragmentation vary with the paleoenvironmental setting? To tackle this question, we compare the size of each morphotype The microfacies ranking corresponds to a ranking of the samples along a theoretical water depth transect, based on microfacies description (Casier et al., 2002;Girard et al., 2014Girard et al., , 2017. Ages are given by the conodont zones (Ziegler and Sandberg, 1990). L.: lower, M.: middle, and U.: upper. The percentages based on counting and cumulated area are provided for the different samples. They correspond to levels in the three outcrops BU, CT, and LS-E . These levels were chosen to document contrasted depositional settings. This was estimated by ranking the samples along a theoretical water depth transect, based on microfacies description in thin sections (Casier et al., 2002;Girard et al., 2014Girard et al., , 2017. To the right: mean area in square micrometers (calculated as the cumulated area of the morphotype divided by the number of item).
to an independent marker of water depth, based on microfacies description of the sediments in thin sections. To do so, we selected a set of samples from three outcrops documenting contrasting water depth.
(2) Do both methods provide a similar result? We compare percentages obtained using the cumulative area to those based on the classical method of counting. Each morphotype was considered separately in a univariate way. Furthermore, the compositions of the assemblages based on area and counting were compared in a multivariate way using a co-inertia analysis.
(3) The compositions of the assemblages based on area and counting are compared to the microfacies description in order to assess if both methods provide a similar relationship to paleoenvironmental estimates.

Material
Samples from three sections were chosen for the study of Late Devonian agglutinated foraminifera (Table 1). Two sections are located in the Montagne Noire (France): the Col des Tribes (CT) and the La Serre trench E (LS-E ) sections. One is located in Thuringia (Germany): the Buschteich (BU) section. Sedimentological studies have shown that CT and BU sections were deposited in an outer carbonated ramp (Girard et al., 2014(Girard et al., , 2017, LS-E being deposited in a more proximal environment (Casier et al., 2002). Along each section, some levels were selected to document contrasted depositional environments. Seven levels were sampled along the CT section and seven along the BU section, and one level was sampled along LS-E . These levels have been attributed to five different facies, based on microfacies description (Girard et al., 2014(Girard et al., , 2017, which have been ranked along a theoretical depth transect (Fig. 2). This provided a rough estimate of the setting in which the foraminifera lived and allowed the testing of the influence of wave actions on breakage.
Agglutinated foraminifera were obtained by dissolving sediments using formic acid (Holcova, 2004), a procedure that precludes documenting calcareous foraminifera and also calcite-cemented agglutinates. This may bias our estimate of the assemblages of agglutinated foraminifers, but in a consistent way in all samples. The term assemblages, used throughout the text, should thus be understood as residual assemblage after formic acid dissolution.
A few grams of rocks were dissolved in 10 % formic acid. After 24 h the residues were washed over a 100 µm sieve. The residues were dried and all foraminifera were manually picked using a stereomicroscope.

Assignment of foraminiferal groups
Based on the outer morphology of their test, foraminifera were attributed to morphotypes. We used the classification of Holcova and Slavik (2013) developed for Middle Devonian foraminifera, complemented for one multilocular morphotype by the classification of Nagy et al. (2009). This corresponds to a simplified version of a morphogroup system developed for more recent foraminifera, from modern to Mesozoic periods (Murray, 2006;Murray et al., 2011;Nagy, 1992). The attribution to morphogroups is based on the general external test morphology. It is therefore independent of the determination at the genus or even species level, based on more detailed features of the test. Although a significant amount of potential ecological information is lost, this technique allows transposition of the interpretation related to morphogroups to other periods of time.
The following five morphotypes were found in the documented samples (Fig. 1). They were defined based on Holcova and Slavik (2013) and Nagy et al. (2009).
-Morphotype ED1 (Fig. 1a, b). It corresponds to tubular foraminifera, either straight (panel a) or irregularly (panel b) meandering. It is interpreted as adapted to an  Table 1) (Casier et al., 2002;Girard et al., 2014Girard et al., , 2017. Deposits are organized according to water depth. BU levels in red, CT levels in green and LS-E levels in black.
erect epifaunal habitat, with a suspension feeding strategy.
-Morphotype ED2. It includes plano-convex forms, which can be hemitubular. It is interpreted as attached epifauna, with a passive herbivore feeding strategy.
-Morphotype ED3. It consists of discoidal coiled tests. The coiling can be irregular of planispiral. It is interpreted as epifaunal, with an active herbivore and/or detritivore feeding strategy.
-Morphotype ED4. It includes globular, spherical forms, which may be simple or with tubular projections. It is interpreted as epifauna to shallow infauna, being a passive deposit feeder.
-Morphotype ED5 . It is not described in Holcova and Slavik (2013), but it corresponds to multilocular forms, as their morphotype ED5. It can be related to subgroups 2b or 4a of Nagy et al. (1992), morphogroup A-4 of Tyszka (1994), or B3 or D of Murray et al. (2011). This morphotype is represented by forms with an initial planispiral coiling, prolongated by a tubular apertural zone. It is interpreted as shallow infauna and an herbivorous and detritivorous feeder.
In each sample, foraminifera were sorted by morphogroup and counted accordingly. The area of each foraminifer was estimated using an image analysis device (Nikon NIS-Elements software) based on pictures taken with the stereomicroscope. Individual areas were cumulated per morphotype. The cumulated area of each morphotype was expressed as the proportion of the total area of the agglutinated foraminifera in a given sample. This provided an area-based estimate of the proportions between morphotypes.

Methods
To assess if the fragmentation was dependent on the depositional context, the area per morphotype (cumulated area of a given morphotype in a given sample divided by the number of items of this morphotype in the sample) was compared to the water depth during deposition, assessed by a ranking of the microfacies along a theoretical gradient of water depth (Casier et al., 2002;Girard et al., 2014Girard et al., , 2017. As this score is a discrete variable, Spearman rank order correlations were used. Note that benthic foraminifers are known to vary with many other aspects of the depositional context, characterizing their microhabitat (e.g., Gooday, 2003), but these aspects were difficult to assess in such ancient sediments. The percentages of each morphotype based on area and counting were first compared one by one, using a linear correlation (Pearson correlation coefficient R). The variation in the composition of the assemblages in the different samples was summarized by a principal component analysis (PCA) on the variance-covariance matrix of the percentages. Two PCAs were thus performed: one on the percentages based on area estimates and one based on the counting. A co-inertia analysis was used in order to evaluate the concordance between the two PCAs in a multivariate way (function coinertia in ade4; Dray and Dufour, 2007). This approach aims at finding orthogonal vectors (i.e., co-inertia axes) maximizing the sum of squared covariances between two datasets (Dodélec and Chessel, 1994;Dray et al., 2003) allowing their projection in a common space.
An RV coefficient further allowed a numerical estimate of the multivariate correlation between the two datasets. This coefficient corresponds to the sum of the squared covariances between two sets of variables, divided by the total amount of variation in the two sets of variables (Escoufier, 1973). The significance of the association was tested by comparing the observed RV coefficient to a distribution obtained by permuting each row separately in each of the two sets of variables (9999 permutations).
The diversity of the assemblages was estimated using the Shannon diversity index based on counting and area (Shannon, 1948). Given p k the frequency of the morphotypes in the assemblages, the Shannon diversity index is calculated as follow: H S = − p k ln(p k ). The two estimates of H S , based on counting and cumulated area, were compared using a linear correlation (Pearson correlation coefficient R). The percentages of the different morphotypes, based on area and counting as well as the Shannon diversity indices, were compared to the microfacies ranking using a Spearman rank order correlation. All statistics were performed using R (R Core Team 2017). The multivariate analyses, including estimates of the RV coefficient and the co-inertia analysis, were performed using the R package ade4 (Dray and Dufour, 2007).

Results
The counting and area estimates (Table 1) provided two images of the composition of the assemblages across the different samples (Fig. 2). Most samples were dominated by the tubular form ED1. The globular ED4 occasionally represented up to 60 % of the assemblage, while the plano-convex ED2, and the multilocular ED5, represented up to 20 % of the assemblages. The discoidal morphotype ED3 appeared to be extremely rare, being only present in one sample in the Buschteich sample BU095. Statistics were therefore not possible for this morphotype, but it was kept in the analyses of the assemblages (percentages, co-inertia, diversity indices).

Variations in the area of the morphotype
The morphotypes varied greatly in their mean size per item, from 20 000 to 160 000 µm 2 (Table 1). None were consistently greater or smaller than the others.

Correlation between compositional estimates based on area and counting
In all cases, the percentages of a morphotype based on area or counting were highly correlated (Pearson product-moment correlation R: ED1 R = 0.9605; ED2 R = 0.9571; ED4 R = 0.9580; ED5 R = 0.9395; in all cases P < 0.001).

Multivariate correlation between compositional estimates based on area and counting
The variations in the assemblages were summarized using PCAs. Considering counting-based percentages, the first axis explained 81.9 % of variance, whereas the second one explained only 13.2 %. The variance explained by the first two axes of the PCA on the area-based percentages were similar (PC1 = 83.0 %, PC2 = 10.5 %). In both cases, the first axis was driven by a balance between ED1 on one side and ED2 and ED4 on the other side. The second axis was driven by ED5 . The contribution of the rare morphotype ED3 is close to the origin, meaning that it plays a small role in structuring the assemblages. In order to compare the ordination of the samples in these two PCAs, we performed a co-inertia analysis (Fig. 3). Coinertia analyses provide a visual representation of how two datasets co-vary, for instance here assemblage compositions based on area or counting. Congruence between the two datasets corresponds to short arrows, and incongruence to long arrows. In the present case, the two datasets match very well, especially regarding the ordination along the first axis, representing 97.7 % of the total projected inertia.
The multivariate correlation between the compositions based on counting and area was further assessed using an RV coefficient. The two datasets were highly correlated (RV = 0.9180, P < 0.001).

Diversity indices based on counting and area
Shannon diversity indices based on counting and area were highly correlated (R = 0.9512, P < 0.001), showing that both methods provide a similar estimate of the assemblage diversity.
Finally, Shannon diversity indices based on counting and cumulated area were not correlated with microfacies score (counting: P = 0.3390; area: P = 0.1990).

Estimating the composition of foraminiferal assemblages: counting vs. cumulated area
Any paleoenvironmental and/or paleoecological interpretation of foraminiferal assemblages first requires reliable estimates of the contribution of the different morphotypes. Tubular forms appear to be one of the key components of some agglutinated foraminiferal assemblages corresponding to deepwater, low-energy habitats (Murray et al., 2011;Setoyama et al., 2011). However, their elongated shape makes them fragile and prone to fragmentation; counting each fragment as a unit may overestimate the contribution of this morphotype to the assemblages (Setoyama et al., 2013). It has thus been proposed to take into account the cumulated length of the tubular fragments to get a better estimate of their abundance (Setoyama et al., 2011). The drawback of this approach is to deliver an estimate for the tubular (ED1) morphogroup only, which is not comparable with other morphotypes estimated by counting, preventing the estimation of the relative proportion of each morphotype. We thus tested here an alternative approach, estimating the cumulated area for each morphotype: this allows an estimate of the composition of the whole assemblage to be compared to the one based on counting. Based on our sampling, documenting a broad range of paleoenvironmental conditions, results based on counting and cumulated area appear to be very close, regarding altogether the proportion of the morphotypes, diversity estimates, and the multivariate composition of the assemblages.
This suggests that counting tubular fragments as units does not significantly overestimate their contribution, a result in agreement with attempts based on considering cumulative length (Setoyama et al., 2011). This further suggests that counting, which is less time-consuming than measuring the cumulated area of each morphotype, may represent a valid method for the estimation of assemblage composition of Late Devonian foraminifera.

Counting and cumulated area: strengths and weaknesses
In most Late Devonian sediments, agglutinated forms are relatively abundant. They can be easily isolated from hardrock sediments using acid attack (Holcova and Slavik, 2013). This method may be responsible for most of the fragmentation and explain the absence of any correlation with the depositional context. If fragmentation was related to taphonomic processes, variations with water depth would have been expected. Increased fragmentation would be expected in the most proximal environments in which surface sediments were subjected to storm wave action, although fragmentation may also increase as well in deeper environments, resulting from low sedimentation rates increasing exposure times. In the case of acid-extracted assemblages, because the bias due to fragmentation is independent of the environment, counting and cumulative area may perform equally well to assess the composition of the assemblages. This remains to be assessed for other ways to isolate benthic foraminifera from the sediment. Cumulative area may be itself prone to bias. An increased cumulated area may not only result from an increased number of specimens in the assemblage but also from specimens of larger size. This may bias the proportions in favor of one morphotype larger than the others. Cumulative area may then provide relevant information about the biomass produced by each morphotype. Intra-morphotype size variations may also occur and could be estimated by dividing the cumulative area by the number of specimens of a given morphotype. Such average size per item may also depend on the environment and thus interfere with the estimation of changes in the assemblage composition (Holcova, 2004). In recent planktic foraminifera, some species grow larger when close to the ecological optimum of the species (Schmidt et al., 2003). Some others grow larger in unfavorable conditions, delaying reproduction and prolonging the ontogenetic life span (Mojtahid et al., 2015;Schiebel and Hemleben, 2017). Such intra-morphotype size variation has been suggested for benthic foraminifera as well (Kaminski and Kuhnt, 1995;Holcova, 2004;Holcova and Slavik, 2013). Assessing cumulative area may be of interest as the variation in the individual size through time and space may be of interest per se in further paleobiological studies, but its estimate will be rendered difficult by the pervasive problem of fragmentation. Using alternate methods for releasing foraminifera from the sediments, using acid but no sieve, may then be of use to minimize fragmentation.

Assemblage composition and microfacies variations
We designed our sampling to maximize the paleoenvironmental variation, in order to compare the relative performance of counting-and area-based methods of assessing the assemblage composition. This sampling, although limited, provided information about the relationship of the morphotypes with the paleoenvironmental setting.
The use of gross morphology to place species into morphogroups that may correspond to ecological preferences has been critically discussed (e.g., Sen Gupta and Machain-Castillo, 1993;Jorissen et al., 2007;Kender and Kaminski, 2017), and the main concerns are that morphogroups do not capture the true ecological preferences of many modern forms and may miss detailed ecological signals. However, we do note that some first-order ecological variations may be captured by morphogroup analysis. Therefore, tubular forms (ED1) are interpreted as epifaunal suspension feeders, favoring deep low-energy environments, with gentle currents sustaining food resource input (Murray et al., 2011). They are frequently interpreted as an indicator for deposition in deep water (e.g., Murray et al., 2011;Nagy et al., 1995;Setoyama et al., 2011) but can also be associated with periods of sealevel fall (Reolid et al., 2012). This complex relationship with water depth is confirmed by the absence of any convincing relationship with the microfacies among our samples. This does not exclude a dependency on the sedimentological and water depth conditions, but it may occur according to a nonlinear relationship that will be difficult to provide evidence for and interpret.
In our sampling, the only morphotype to display a relationship to the depositional context was the multilocular form ED5 , interpreted as shallow infauna, herbivorous to active detritivorous scavenging (Nagy et al., 2009). In modern environments, this form seems to be dominant in relatively proximal environments of shelf (Jones and Charnock, 1985;Murray, 2006), an interpretation that would fit their increase in relative abundance in the most proximal conditions. Ecological preferences of the three other morphotypes remain to be better assessed.
This challenges future studies based on more extensive sampling through time and space to better understand the environmental controls on the composition of benthic foraminifera in the Late Devonian period and other time periods and environmental settings.

Conclusions
Our results show that both methods of estimating morphotype percentages, the traditional counting and the cumulated area methods, provide similar results, are highly correlated with each other, and provide similar relationships with paleoenvironmental proxies. We find no clear relationship between the depositional setting and fragmentation, which suggests that most breakage occurs during the dissolution of rocks. The induced bias should thus be similar across samples. Overall, this suggests that the classical method based on counts, which is less time-consuming than measuring the cumulated area of each morphotype, is only partially biased by the issue of fragmentation. This implies that traditional counting methods are an appropriate way of capturing foraminiferal assemblage characteristics, even when specimens are fragmented. Data availability. All data are available in the Table 1.