Philippe Janssen
email: philippe.janssen@irstea.fr

Stéphane Bec
email: stephane.bec@univ-grenoble-alpes.fr

Marc Fuhr
email: marc.fuhr@irstea.fr

Pierre Taberlet
email: pierre.taberlet@univ-grenoble-alpes.fr

Jean-Jacques J Brun
email: jean-jacques.brun@irstea.fr

Christophe Bouget
email: christophe.bouget@irstea.fr

Present conditions may mediate the legacy effect of past land-use changes on species richness and composition of above-and below-ground assemblages

Keywords: ancient forest, community dynamics, dispersal, forest management, habitat quality, land-use history, mountain forest, plant-soil interactions, secondary succession

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The influence of landscape history on contemporary biodiversity patterns is increasingly recognized as a key factor that structures species assemblages today [START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF]. Given that most of the world's terrestrial ecosystems have been more or less impacted by human actions [START_REF] Newbold | Global effects of land use on local terrestrial biodiversity[END_REF], understanding legacies of the past may help to better understand present biodiversity patterns and predict the future ecological impacts of on-going human practices on ecosystem services and functions [START_REF] Vellend | Historical ecology: Using unconventional data sources to test for effects of global environmental change[END_REF]. In a context of continuing global biodiversity loss [START_REF] Butchart | Global Biodiversity: Indicators of Recent Declines[END_REF], understanding the influence of past and present human actions on ecosystems is of primary importance for conservation.

Forest ecosystems represent one third of Europe's total land area (Forest Europe 2015). However, European forest cover has fluctuated widely over the last centuries [START_REF] Kaplan | The prehistoric and preindustrial deforestation of Europe[END_REF].

After a long period of deforestation, forest cover has consistently increased since the mid-nineteenth century (FAO 2015) and has nearly doubled in France [START_REF] Mather | The course and drivers of the forest transition: the case of France[END_REF]. This reforestation has mainly occurred on abandoned agricultural land, either through spontaneous recolonization or deliberate replanting.

Two types of forests have resulted: ancient forests, which have existed continuously for centuries, and recent forests, which result from reforestation after a certain threshold date [START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF]. This massive land-use change has led researchers to develop the concept of "forest continuity", also referred as "forest ancientness" or "ancient woodland" [START_REF] Flinn | Recovery of forest plant communities in post-agricultural landscapes[END_REF][START_REF] Goldberg | The ancient woodland concept as a practical conservation tool in Great Britain[END_REF][START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF]).

Forest continuity is thus defined as a minimum residence time of the wooded state since a threshold date, which differs between countries in northwestern Europe due to the complex land use history and availability of historical maps (e.g. 1600 or 1700 in GB; 1750 or 1800 in Germany; 1770-1800 in Belgium & Denmark;1820 in Sweden;1820-1850 in France & Netherlands;[START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF].

Forest continuity has been shown to play a fundamental role in driving herbaceous plant assemblages in temperate forests [START_REF] Flinn | Recovery of forest plant communities in post-agricultural landscapes[END_REF][START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF]. Two processes have been highlighted: dispersal limitations due to poor ability of ancient-forest plant species to colonize recent forests (e.g. [START_REF] Verheyen | An integrated analysis of the effects of past land use on forest herb colonization at the landscape scale[END_REF][START_REF] Naaf | Colonization credit of post-agricultural forest patches in NE Germany remains 130-230years after reforestation[END_REF] and recruitment limitations due to modifications in soil properties and competitive interactions (e.g. Honnay et al. 2002a;[START_REF] Baeten | Environmental limitation contributes to the differential colonization capacity of two forest herbs[END_REF]. Given the influence of forest continuity on primary producers it is therefore reasonable to assume that continuity has an indirect effect on primary consumers and decomposers, especially for the less motile taxa. For example, it has been shown that changes in plant communities between ancient and recent forests act indirectly on true bugs [START_REF] Gossner | Plant and arthropod communities in young oak stands: are they determined by site history?[END_REF]. Moreover, due to modifications in soil properties related to past land-use, forest continuity can be expected to have direct effects on below-ground soil fauna. However, very few studies to date have investigated these relationships (e.g. [START_REF] Assmann | The ground beetle fauna of ancient and recent woodlands in the lowlands of north-west Germany (Coleoptera, Carabidae)[END_REF][START_REF] Buse | Ghosts of the past'": flightless saproxylic weevils (Coleoptera: Curculionidae) are relict species in ancient woodlands[END_REF][START_REF] Janssen | Are biodiversity patterns of saproxylic beetles shaped by habitat limitation or dispersal limitation? A case study in unfragmented montane forests[END_REF].

Beyond past land-use effects, present human activities also impact forest biodiversity. Through management practices, foresters substantially change stand ecological conditions, and this has consequences on biodiversity [START_REF] Paillet | Biodiversity Differences between Managed and Unmanaged Forests: Meta-Analysis of Species Richness in Europe[END_REF]. Indeed, regardless of forest continuity, many forest taxa are influenced by stand maturity attributes, e.g. dead-wood volume. This suggests that there is an additive effect of stand maturity over forest continuity which, if not accounted for, may mask or reinforce the expected ecological processes studied. However, very few studies have controlled for stand maturity when comparing ancient and recent forests (but see, [START_REF] Kolb | Effects of environment, habitat configuration and forest continuity on the distribution of forest plant species[END_REF][START_REF] Kelemen | Effects of land-use history and current management on ancient woodland herbs in Western Hungary (ed B Collins)[END_REF], making it difficult to distinguish the effect of stand maturity from the effect of forest continuity per se [START_REF] Nordén | Effects of Ecological Continuity on Species Richness and Composition in Forests and Woodlands: A Review[END_REF]). To properly disentangle the effect of past and present human activities and to subsequently improve conservation strategies, this information is fundamental.

We aimed to study the effect of forest continuity on vascular plants and springtails in mixed montane forests, while controlling for stand maturity. Since forest continuity is expected to influence plant assemblages and since plant assemblages determine resource inputs for the decomposer subsystem [START_REF] Hansen | Effects of habitat complexity and composition on a diverse litter microarthropod assemblage[END_REF][START_REF] Wardle | The influence of biotic interactions on soil biodiversity[END_REF]), variations in plant assemblages between ancient and recent forests may have important effects on below-ground biota and the processes that they regulate [START_REF] Wardle | Ecological Linkages Between Aboveground and Belowground Biota[END_REF]. Given the huge number of diverse species and the fundamental role played by soil fauna in ecosystem functioning [START_REF] Decaëns | The values of soil animals for conservation biology[END_REF], it appear essential to investigate the legacy effects induced by forest continuity on below-ground biodiversity.

Moreover, due to their small body size and the difficulties they have moving within the soil, most soil taxa are likely to have poor dispersal abilities. As such, it has been showed that springtail diversity is affected by landuse change and that species recovery in reforested sites may take more than 50 years, especially for poorly dispersing species (Ponge et al. 2006). Moreover, springtails respond to modifications in soil conditions and vegetation cover [START_REF] Ponge | Collembolan communities as bioindicators of land use intensification[END_REF][START_REF] Korboulewsky | How tree diversity affects soil fauna diversity: A review[END_REF], two ecosystem alterations that may result from past and present human-induced activities. Thus, a combined above-and below-ground approach should enhance our understanding of the human impacts on biodiversity, and result in a better orientation of conservation strategies [START_REF] Wardle | Ecological Linkages Between Aboveground and Belowground Biota[END_REF]. Specifically, we tested whether species richness and composition of herbaceous plants and springtails responded to forest continuity and stand maturity. As it is likely that shifts in species richness and composition result from a complex of different factors acting at different spatial scales [START_REF] Siefert | Scale dependence of vegetation-environment relationships: a meta-analysis of multivariate data (ed A Acosta)[END_REF], we also tested the relative influence of a set of soil, stand-and landscape-scale variables on plants and springtails. Landscape variables were used to consider possible factors underlying the effect of forest continuity, e.g. related to differences in the amount of habitats in the surrounding landscapes between ancient and recent forest sites.

Stand and soil variables were used to account for important local characteristics that may structure plant and springtail communities, beyond the influence of stand maturity and forest continuity. Moreover, soil variables helped us to control for possible remaining effects of past land-uses (e.g. [START_REF] Dupouey | Irreversible impact of past land use on forest soils and biodiversity[END_REF][START_REF] Dambrine | Present forest biodiversity pattern in France related to former roman agriculture[END_REF]). Based on this scheme, we addressed the following two questions: (i) Are past human-induced activities or are present human-induced activities more important in shaping herbaceous plant and springtail species richness and composition? (ii) Among soil, stand and landscape, which factor has the greatest influence on herbaceous plant and springtail species richness and composition?

Materials and methods

Study area and experimental design

The study was carried out in the French pre-Alps in the Vercors, Chartreuse and Bauges ranges (Fig. 1), all of which are characterized by a limestone substratum and a temperate climate. The landscapes are mostly forested (63%) and reforestation has mainly occurred above and below the persistent forest belt, as it has elsewhere in European mountain areas [START_REF] Kozak | Forest Cover Change in the Western Carpathians in the Past 180 Years: A Case Study in the Orawa Region in Poland[END_REF][START_REF] Gellrich | Agricultural land abandonment and natural forest re-growth in the Swiss mountains: A spatially explicit economic analysis[END_REF]. Moreover, due to physical constraints and lack of logging roads, mountain forests has hitherto been less intensively managed than lowland forests [START_REF] Paillet | Quantifying the recovery of old-growth attributes in forest reserves: A first reference for France[END_REF]. Therefore, compared to recent lowland forests, the recent montane forests in the Northern Alps, i.e. forests that have existed since at least 1864, are mostly adjacent to ancient forests and have the potential to develop towards stand structures similar to those found in ancient forests.

We sampled 70 sites located in montane beech-fir forests at an altitude of 800 -1500 m. The dominant tree species are European beech Fagus sylvatica, silver fir Abies alba and Norway spruce Picea abies.

Our stratified sampling design crossed forest continuity (ancient forests = 37; recent forests = 33) and stand maturity (mature stands = 37; overmature stands = 33), resulting in 22 ancient-mature sites, 15 ancientovermature sites, 15 recent-mature sites and 18 recent-overmature sites. All selected sites were managed by selective cutting at some point, i.e. stumps were measured at and around sites (see Appendix S3 & S4). To ensure the independence among observations and avoid edge effects, all the sampling sites were established > 1 km away from any other sites, were located in ancient or recent forests > 5 ha in area, and were > 50 m from the nearest stand edge.

Forest continuity, stand maturity and soil characterization

Forest continuity was characterized by crossing digitized and geo-referenced 1:40 000 État-Major maps of France, charted in the middle of the 19 th century, combined with up-to-date 1:10 000 vegetation maps in a Geographic Information System managed with QGIS (QGIS Development Team 2015). Forest cover overlapping in both maps was considered to indicate ancient forests, while current forest cover overlapping with crops, pastures or meadows in the État-Major maps was considered to indicate recent forests. Around each selected ancient and recent forest site (500 m radius), we controlled the État-Major maps using 1:5 000 or 1:2 500 cadastral plans drawn between 1809 and 1838 for the Vercors and the Chartreuse ranges ("Napoléonic cadastre") and between 1730 and 1738 for the Bauges range ("Sarde cadastre"). Past land-use in recent forest sites was then accurately described (pastures: n = 17; meadows: n = 13; crops: n = 3). Moreover, to confirm the continuity of the forest cover since the middle of the 19 th century, we used aerial photographs taken in the 1950s, at and around each site.

Stand maturity was a priori approached by forest prospections (n = 117) and rapid stand attribute measurements (diameter of the three largest trees, presence of large snags and/or of large coarse woody debris). We then selected 70 stands in a way that ensures a balanced distribution of sites between ancient and recent forests and a priori mature and overmature stands. In each selected stands, stand maturity was then characterized on a 20 m radius plot and a 10 m radius subplot in which all standing trees and lying trunks were recorded (for further details, see [START_REF] Janssen | Are biodiversity patterns of saproxylic beetles shaped by habitat limitation or dispersal limitation? A case study in unfragmented montane forests[END_REF]. Canopy openness was estimated with a spherical densiometer from four points in the cardinal directions, 10 m away from the plot center. To distinguish between mature and overmature stands (see Appendix S2), we used hierarchical cluster analysis (Ward method) based on the three first axes of a principal component analysis (PCA) with a cumulative projected inertia of 88.30% (ade4 package). PCA was conducted considering four environmental variables closely related to stand maturity: volume of large coarse woody debris (Ø > 30 cm), number of large snags (Ø > 30 cm), number of very large living trees (Ø > 62.5 cm) and tree microhabitat diversity (i.e. the number of microhabitat types: cavities, sporophores of saproxylic fungi, ivy, sap runs, missing bark, cracks and shelter bark).

To characterize soil properties, we dug a soil pit > 50 cm in depth in the plot center. In each soil pit, two soil samples were collected -at 0-10 cm (topsoil) and at 30-40 cm (subsoil) -and humus forms were described. Each soil sample was analysed for pH (H 2 O, NF ISO 10390), total C (NF ISO 10694), total N (NF ISO 13878) and available P content (Olsen method, NF ISO 11263). Subsoil samples were also analysed for clay, silt and sand proportions. Soil analyses were performed by the INRA Soil Analysis Laboratory following procedures of the French Standardization Association (i.e. AFNOR).

Environmental variables

To model the response of plants and springtails to forest continuity and stand maturity, we used soil, stand and landscape variables in addition to these two factors. Among a larger set of potential variables (see (H_index), the pH (pH) and the total carbon/total nitrogen ratio (CN_ratio) -were derived from humus form description and topsoil sample analysis.

Appendix

Plant survey

Following the Braun-Blanquet abundance-dominance methodology, we surveyed all the vascular plants within the 10 m radius circular subplot. To limit variation in plant species composition due to environmental conditions, all the plots were established in ecologically similar sites (mesophilic and neutrophilic) belonging to the alliance Fagion sylvaticae. Moreover, in order to characterize the entire groundlayer plant community and account for plant phenology, two complete vegetation surveys were done: from May to June 2014 for vernal species and from July to August 2014 for estival species. Finally, to limit biases associated with observer effect, all surveys were performed by the first author (PJ).

Springtail assessment

Soil fauna studies often suffer from the inefficiency of sampling strategies and a strong taxonomic deficit [START_REF] André | Soil biodiversity: myth, reality or conning?[END_REF]. We therefore chose to use environmental DNA (eDNA) from the soil core samples to characterize springtail communities. eDNA is an emerging approach for biodiversity studies [START_REF] Thomsen | Environmental DNA -An emerging tool in conservation for monitoring past and present biodiversity[END_REF] that has already proved its efficiency, e.g. for earthworms [START_REF] Pansu | Landscape-scale distribution patterns of earthworms inferred from soil DNA[END_REF].

In each site, two semi-circular 10 m radius plots were established. Within each semi-circular plot, 25 soil cores (10-cm depth) were collected in September 2015 and pooled together. To avoid cross-contamination, the coring sampler was sterilized between sites by flame cleaning and, to avoid DNA degradation before extraction, soil samples were dehydrated with silica gel. For each composite soil sample (n = 140), soil extracellular DNA was isolated in November 2015 following Taberlet et al.'s (2012) protocol. Extraction was performed with the NucleoSpin® Soil isolation kit (Macherey-Nagel, Düren, Germany); the lysis step was skipped and the manufacturer's instructions were followed closely.

Springtail communities were assessed with a primer pair targeting a region of the 16S rRNA gene. For each sample, three PCR replicates were carried out, resulting in 420 samples (see Appendix S1). Library construction and sequencing for the samples (Illumina HiSeq 2500 platform, 2 x 100 bp, paired-end reads) were carried out at Fasteris (Geneve, Switzerland). The DNA sequences read were then filtered and clustered into molecular operational taxonomic units (MOTU) with the dedicated OBITools package and according to an established workflow [START_REF] Boyer | OBITOOLS: a UNIX-inspired software package for DNA metabarcoding[END_REF]. A complete description of data processing is proposed in Appendix S1.

Statistical analysis

Analyses were performed with R version 3.2.2 (R Core Team 2016). We used presence/absence for plants and for springtails as dependent variables in the following analysis. Based upon data exploration (Zuur, Ieno & Elphick 2010), independent variables with a skewness >1 were log or log+1 transformed to approximate normal distributions. For proportional data, logit transformation was applied [START_REF] Warton | The arcsine is asinine: the analysis of proportions in ecology[END_REF]. We then used two-way ANOVAs with type III sum of squares to test the significance of each environmental variable to forest continuity, stand maturity and their interaction term (Table 1, see Appendix S3 & S4).

To determine whether plant or springtail richness was influenced by forest continuity and/or stand maturity, we also used two-way ANOVAs with type III sum of squares. To verify hypothesis statements, we developed 49 a priori candidate GLMs, plus a null model (see Appendix S6). We used Poisson regressions to relate richness in plants and springtails to soil, stand and landscape variables; we considered a balanced number of models to avoid biasing estimates toward a particular category of environmental variables. We assessed multi-colinearity among explanatory variables with variance inflation factors (car package) and spatial autocorrelation with Moran's correlogram in the model residuals among sampling sites (ncf package). To identify the most parsimonious GLM, we used Akaike's information criterion corrected for small sample sizes (AICc, [START_REF] Burnham | Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach[END_REF]. Moreover, we used model averaging to estimate parameter and associated unconditional standard errors based on the subset of top ranking models for which the sum of the AICc weights reached ≥ 0.95 (MuMIn package). To determine the relative importance of soil versus stand models, soil versus landscape models and stand versus landscape models, we compared the sum of the AICc weights for three categories of models: e.g. models that accounted for soil variables only, for stand variables only and models that accounted for both types of variables [START_REF] Burnham | Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach[END_REF].

To determine if species composition was influenced by forest continuity and/or stand maturity, we used multivariate GLMs. This approach has the advantage of being based on a GLM framework, which is typically better able to handle the strong mean-variance relationship in abundance data [START_REF] Warton | Distance-based multivariate analyses confound location and dispersion effects: Mean-variance confounding in multivariate analysis[END_REF]. We fitted separate models to each species with a binomial distribution and summed across the univariate response to estimate their multivariate response, i.e. assuming the independence of the species response variables. The significance of the multivariate response was assessed with the PIT-trap method with 999 bootstrap resamples (mvabund package). We first assessed the effect of forest continuity and/or stand maturity on plant and springtail composition and then performed a canonical analysis of principal coordinates (CAP, [START_REF] Anderson | Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology[END_REF] with a Sørensen-Dice index, to provide a graphical representation of the interaction. We then ran the 49 a priori models, plus the null model (see Appendix S6). To identify the most parsimonious multivariate GLM, we used AICc. For each multivariate model, AICc was obtained by extracting the log-likelihood and calculating the AICc of each univariate model, summing AICc values and standardizing the sum by the number of observations. Model performance was evaluated according to the proportion of deviance explained (D 2 = (null deviance -residual deviance) / null deviance), and the significance of each independent variable was assessed using an analysis of variance on the full multivariate GLM. To determine the relative importance of soil versus stand models, soil versus landscape models and stand versus landscape models, we compared the sum of the AICc weights for the three categories of multivariate GLMs. ). Indeed, nearly all of the variables related to stand structural complexity and dead-wood compartment increased from mature to overmature stands; while all the variables related to forest proportion and distance to the edge increased from recent to ancient forests. On the other hand, soil variables

Results

Overall

were not clearly related to forest continuity or to stand maturity. Only some weakly significant relationships were found: topsoil C/N ratio and subsoil clay proportion increased from ancient to recent forests, while total subsoil nitrogen increased from recent to ancient forests; humus index and total subsoil carbon increased from overmature to mature stands; while subsoil clay proportion increased from mature to overmature stands.

However, total phosphorous and total nitrogen did not increase and pH did not become more basic in either subsoil or topsoil when ancient and recent forests were compared, indicating the absence of a legacy effect related to past land-use.

The interaction term between forest continuity and stand maturity was non-significant for almost all the tested variables (see Appendix S4), indicating that these variables varied consistently between ancient and recent forests, at a comparable level of maturity.

Response of species richness to forest continuity, stand maturity and associated variables

Two-way ANOVAs revealed no difference in plant or springtail richness between ancient and recent forests, mature and overmature stands or between ancient-mature, ancient-overmature, recent-mature and recent-overmature stands.

GLM results showed that plant and springtail richness were best predicted by the same most parsimonious model (i.e. based on AICc), which accounted for soil and stand variables (Table 2). Observed

Moran's I for the top model's residuals was non-significant for both dependent variables, indicating that spatial patterns had been accounted for by the environmental variables (see Appendix S7). The predictive value of this model was high for plants (pseudo-R² = 73%) and low for springtails (pseudo-R² = 16%). Although this model seems quite robust for explaining richness patterns, model selection uncertainty still remains since associated AIC c weight and evidence ratio were relatively low. We therefore used model averaging -the 7 top-ranked plant models and 38 top-ranked springtail models -to draw inferences about the variables' influence on richness. Plant and springtail richness increased with an increasing pH value (Table 3). Stand variables only influenced plant richness, which increased with mean canopy openness and the proportion of conifers but decreased with increasing number of tree diameter classes. Landscape variables did not influence either plant or springtail richness.

Judging from the sum of AICc weights for GLMs, plant richness was more influenced by models accounting for a combination of different environmental factors than was springtail richness (Fig. 2). Indeed, for soil versus stand models, plant richness was better explained by models combining both types of variables (sum of AIC c weight, ∑AICw = 0.87), while springtail richness was more heterogeneous, though somewhat better explained by soil models (∑AICw = 0.44). When comparing soil versus landscape models, both plant and springtail richness was better explained by soil models only (∑AICw, range 0.56 -0.67). Finally, when comparing stand versus landscape models, plant richness was slightly better explained by models combining both types of variables (∑AICw = 0.50), while springtail richness was more heterogeneous, though better explained by stand models (∑AICw = 0.41).

Response of species composition to forest continuity, stand maturity and associated variables

Multivariate GLMs showed that the interaction term between forest continuity and stand maturity was significant for both plant (D² = 0.099, p = 0.001) and springtail (D² = 0.076, p = 0.021) assemblages.

Variation in species composition between ancient and recent forests were thus different in mature and overmature stands. Constrained ordination plots showed that plant and springtail species composition converged from mature to overmature stands, regardless of forest continuity (Fig. 3). In others words, dissimilarities among assemblages of species were more important between ancient and recent forests in mature stands.

Multivariate GLMs results showed that variations in plant and springtail species composition were best predicted by the same most parsimonious model (i.e. based on AICc), which accounted for soil variables only (Table 4). However, compared to the full model (plants, D 2 = 40%; springtails, D 2 = 28%), this model explained a relative low proportion of the variance. Moreover, model selection uncertainty remains for both dependent variables since associated AIC c weight and evidence ratio were low and since the null model was one of the ten top-ranked models. Regarding the three most important variables influencing composition patterns (based on the proportion of deviance explained), our results showed that plants were especially influenced by humus index, the number of tree diameter classes and the proportion of coniferous; while springtails were especially influenced by humus index, pH and canopy openness (Table 5).

Judging from the sum of AICc weights for multivariate GLMs, plant and springtail species composition was more heterogeneous and models accounting for a combination of different environmental factors were of relatively low importance (Fig. 4). Indeed, for soil versus stand models, plant assemblages were equally well explained by both soil and stand models, while springtail assemblages were better explained by soil models only (∑AICw = 0.54). When comparing soil versus landscape models, both plant and springtail assemblages (∑AICw= 0.57 and ∑AICw = 0.74, respectively) were better explained by soil models. Finally, when comparing stand versus landscape models, plants assemblages were equally well explained by stand and landscape models, while springtail assemblages were better explained by stand models only (∑AICw = 0.68).

Discussion

Our results clearly show that plant and springtail species richness and composition were more influenced by present human-induced activities than by past human-induced activities. Moreover, no legacy effect of past land-use were found for either soil or stand properties. Specifically, we highlighted the influence of site-specific characteristics (i.e. soil properties) and present management practices (i.e. stand structure and composition) in shaping above-and below-ground community structure. Landscape scale factors were not important in explaining variations in plant and springtail species richness and composition. Furthermore, we demonstrated that, in unfragmented montane forests dominated by ancient forests, the effect of forest continuity is mediated by stand maturity and soil properties.

Richness patterns are obviously shaped by local factors

Species richness was not directly influenced by stand maturity, let alone by forest continuity. This may indicate the absence of past and present human-induced activities on current total above-and below-ground richness. Previous studies also reported that forest continuity had no influence on total plant richness (e.g. [START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF][START_REF] Sciama | Floristic and ecological differences between recent and ancient forests growing on non-acidic soils[END_REF]. For springtails, only one study has reported that richness increases in ancient forests, as compared to very recently established forests [START_REF] Heiniger | Effect of habitat spatiotemporal structure on collembolan diversity[END_REF]. For stand maturity, many studies have examined the influence of stand age on plant richness or have compared the plant richness of managed versus unmanaged stands (as a proxy for stand maturity), but no clear trends have been reported [START_REF] Paillet | Biodiversity Differences between Managed and Unmanaged Forests: Meta-Analysis of Species Richness in Europe[END_REF][START_REF] Duguid | A meta-analysis of the effect of forest management for timber on understory plant species diversity in temperate forests[END_REF]. Likewise, for springtails, while some studies state that richness increases with stand age (e.g. [START_REF] Salamon | The Collembola community of pure and mixed stands of beech (Fagus sylvatica) and spruce (Picea abies) of different age[END_REF]Chauvat et al. 2011b), others report a non-significant increase or even a decrease (e.g. [START_REF] Addison | Abundance, species diversity, and community structure of Collembola in successional coastal temperate forests on Vancouver Island, Canada[END_REF][START_REF] Chauvat | Successional changes of Collembola and soil microbiota during forest rotation[END_REF]. Overall our results are consistent with current knowledge and indicate that the classification used, i.e. ancient versus recent and mature versus overmature, may be too coarse and heterogeneous to properly reveal the underlying ecological processes involved. To disentangle the effect of forest continuity and stand maturity on species richness, it therefore appears essential to move beyond simple categories to account for related factors.

As compared to soil and stand scale factors, landscape scale factor performed poorly in explaining richness patterns. This indicates that, at the spatial scale of our study, site-specific characteristics are of great importance in shaping above-and below-ground species richness. However, this does not rule out the possibility that richness patterns are better explained by other landscape factors at larger spatial scales.

Specifically, soil factors were powerful predictors of plant and springtail richness, while stand factors were powerful predictors of plant richness only. Several studies have pointed out the importance of local-scale factors for species richness [START_REF] Wang | Relative importance of climate vs local factors in shaping the regional patterns of forest plant richness across northeast China[END_REF][START_REF] Zellweger | Disentangling the effects of climate, topography, soil and vegetation on stand-scale species richness in temperate forests[END_REF][START_REF] Bernhardt-Römermann | Drivers of temporal changes in temperate forest plant diversity vary across spatial scales[END_REF]. Light availability is known to influence plant richness [START_REF] Paillet | Biodiversity Differences between Managed and Unmanaged Forests: Meta-Analysis of Species Richness in Europe[END_REF][START_REF] Bernhardt-Römermann | Drivers of temporal changes in temperate forest plant diversity vary across spatial scales[END_REF], since semi-open stands allow light-and shade-tolerant species to coexist. Moreover, through variations in closely related resources, i.e. light, water and soil nutrients, overstory structure and composition influence understory richness [START_REF] Barbier | Influence of tree species on understory vegetation diversity and mechanisms involved-A critical review for temperate and boreal forests[END_REF]. For example, the positive influence of the coniferous proportion we found on plant richness may be related to an increase in soil water availability, given that the shift in composition was driven by silver fir, a tree that prefers meso-hygrophilous soil conditions. For springtails, it has been suggested that tree mixture, through the diversity of food resources supplied, may positively influence species richness [START_REF] Salamon | The Collembola community of a Central European forest: Influence of tree species composition[END_REF]Chauvat et al. 2011a). As previously stated [START_REF] Korboulewsky | How tree diversity affects soil fauna diversity: A review[END_REF], our results indicate that stand factors are of secondary importance and that springtail richness is predominantly influenced by soil factors, especially pH [START_REF] Salmon | Changes in humus forms and soil animal communities in two developmental phases of Norway spruce on an acidic substrate[END_REF][START_REF] Salamon | The Collembola community of pure and mixed stands of beech (Fagus sylvatica) and spruce (Picea abies) of different age[END_REF]. Apart from forest-continuity and stand-maturity associated factors, it appears essential to account for soil factors to better predict species richness. Our findings point out the importance of local resource availability and heterogeneity in driving richness patterns of above-and below-ground taxa [START_REF] Bartels | Is understory plant species diversity driven by resource quantity or resource heterogeneity?[END_REF][START_REF] Decaëns | Macroecological patterns in soil communities: Soil community macroecology[END_REF].

Species composition: the effect of forest continuity is mediated by stand maturity and local factors

Species composition was not directly influenced by forest continuity. Indeed, plant and springtail assemblages were closely related in both ancient and recent overmature stands, but strongly differed in ancient and recent mature stands. This indicates that the effect of forest continuity was mediated by stand maturity and thus that present human-induced activities may be of more importance than past ones in shaping assemblages patterns. This surprisingly slight legacy effect is interesting given the large number of studies that have demonstrated the influence of forest continuity on plant assemblages (review in [START_REF] Flinn | Recovery of forest plant communities in post-agricultural landscapes[END_REF][START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF]. However, most of the previous studies were conducted in rather similar ecological conditions: lowland temperate forests encompassing fragmented and forest-poor landscapes and dominated by recent forests [START_REF] De Frenne | Interregional variation in the floristic recovery of post-agricultural forests: Recovery of post-agricultural forest[END_REF]). In the Northern Alps where our study plots were located, the landscapes are mainly covered by unfragmented ancient forests, and reforestation has largely occurred next to these ancient forests. Since fragmentation worsens dispersal limitations [START_REF] Jamoneau | Fragmentation alters beta-diversity patterns of habitat specialists within forest metacommunities[END_REF] and since colonization efficiency depends on habitat availability [START_REF] De Frenne | Interregional variation in the floristic recovery of post-agricultural forests: Recovery of post-agricultural forest[END_REF], forest-plant species have been able to rapidly colonize recent forests in our study area. For springtails, [START_REF] Huhta | Collembolan communities in deciduous forests of different origin in Finland[END_REF] reported that assemblages differed between 30-years-old birch plantations established on spruce forest soils and on arable land. However, the authors found that soil properties differed with forest continuity and were key factors in explaining differences in species composition. In our case, and contrary to numerous studies (e.g. [START_REF] Dupouey | Irreversible impact of past land use on forest soils and biodiversity[END_REF][START_REF] Dambrine | Present forest biodiversity pattern in France related to former roman agriculture[END_REF], no legacy effect was found for soil properties. We therefore considered soil habitat quality to be equivalent between the ancient and recent forests. Moreover, most of our recent forest sites were established on former pastures or meadows. Since the legacy effect of past land-use is more important on former croplands than on former grasslands [START_REF] Koerner | Influence of past land use on the vegetation and soils of present day forest in the Vosges mountains, France[END_REF][START_REF] Wulf | Plant species richness of afforestations with different former use and habitat continuity[END_REF] and since remnant population of forest species are better able to survive in grasslands [START_REF] Wulf | Plant species richness of afforestations with different former use and habitat continuity[END_REF][START_REF] Burst | Interactive effects of land-use change and distance-to-edge on the distribution of species in plant communities at the forest-grassland interface[END_REF], recolonization of recent forests have probably been accelerated. This is especially true in mountain areas where, as compared to lowland, a larger range of herbaceous forest species has the ability to grow in open habitats [START_REF] Schmidt | Waldartenlisten Der Farn-Und Blütenpflanzen, Moose Und Flechten Deutschlands, BfN-Skripten[END_REF]. We then inferred that the plant and springtail assemblages in our study were not influenced by recruitment and dispersal limitations between ancient and recent forests and that, therefore, forest continuity was of minor importance in shaping assemblages patterns. Instead, our results underline the importance of stand maturity. Changes in stand characteristics with stand age or due to forestry practices are well-known to influence both plant (e.g. [START_REF] Scheller | Understory species patterns and diversity in old-growth and managed northern hardwood forest[END_REF][START_REF] Amici | Influence of secondary forest succession on plant diversity patterns in a Mediterranean landscape[END_REF][START_REF] Sabatini | Ground Layer Plant Species Turnover and Beta Diversity in Southern-European Old-Growth Forests (ed B Bond-Lamberty)[END_REF]) and springtail communities (e.g. [START_REF] Addison | Abundance, species diversity, and community structure of Collembola in successional coastal temperate forests on Vancouver Island, Canada[END_REF][START_REF] Chauvat | Successional changes of Collembola and soil microbiota during forest rotation[END_REF][START_REF] Salamon | The Collembola community of pure and mixed stands of beech (Fagus sylvatica) and spruce (Picea abies) of different age[END_REF]. Here, we document for the first time that stand maturity may mediate the effect of forest continuity by leading to a convergence in assemblage patterns with the increasing availability of overmature stand attributes. In recent forests, the effect of stand maturity may also reflect an increase in the residence time of the wooded state from mature to overmature stands. This may have led to an increase of the recolonization potential of forest species in recent-overmature stands but also to a greater recovery in environmental conditions between ancient-and recent-overmature stands (e.g. [START_REF] Flinn | Agricultural legacies in forest environments: tree communities, soil properties, and light availability[END_REF]. Since stand maturity may cause strong environmental filtering on biodiversity, and since stand maturity may interact with forest continuity, controlling for habitat quality seems essential to demonstrate the effect of habitat continuity per se [START_REF] Nordén | Effects of Ecological Continuity on Species Richness and Composition in Forests and Woodlands: A Review[END_REF]; especially in geographical areas where unfragmented forests dominate and where soil properties have little or no legacy effect.

Compared to soil and stand scale factors, landscape scale factors performed poorly in explaining plant and springtail species composition. This confirms the importance of site-specific characteristics in shaping above-and below-ground communities. Although several studies have shown that ancient-forest species distribution is influenced by landscape configuration and composition (Honnay et al. 2002b;[START_REF] De Frenne | Interregional variation in the floristic recovery of post-agricultural forests: Recovery of post-agricultural forest[END_REF][START_REF] Naaf | Colonization credit of post-agricultural forest patches in NE Germany remains 130-230years after reforestation[END_REF], it has also been reported that, in unfragmented forests dominated by ancient forests, local factors explain the larger part of species composition [START_REF] Jamoneau | Fragmentation alters beta-diversity patterns of habitat specialists within forest metacommunities[END_REF]). In the Northern-Alps, species assemblages were indeed largely influenced by local factors. These findings are consistent with previous studies showing the importance of soil pH and humus forms, stand composition and structure or canopy openness in shaping both plant and springtail assemblages (e.g. [START_REF] Scheller | Understory species patterns and diversity in old-growth and managed northern hardwood forest[END_REF][START_REF] Salmon | Changes in humus forms and soil animal communities in two developmental phases of Norway spruce on an acidic substrate[END_REF][START_REF] Salmon | Relationships between soil fauna communities and humus forms: Response to forest dynamics and solar radiation[END_REF][START_REF] Sabatini | Ground Layer Plant Species Turnover and Beta Diversity in Southern-European Old-Growth Forests (ed B Bond-Lamberty)[END_REF][START_REF] Korboulewsky | How tree diversity affects soil fauna diversity: A review[END_REF]. The importance of these local factors points out the necessity to account for site specific characteristics and current human-induced activities when analyzing the effects of forest continuity. Though it has been shown that tree species influences the recruitment of ancient-forest species [START_REF] Thomaes | Tree species determine the colonisation success of forest herbs in post-agricultural forests: Results from a 9 yr introduction experiment[END_REF], that changes in canopy structure and composition induce a shift in plant communities in ancient deciduous forests [START_REF] Verheyen | Driving factors behind the eutrophication signal in understorey plant communities of deciduous temperate forests[END_REF], and that several ancient-forest plant species decrease in abundance with local disturbances [START_REF] Naaf | Traits of winner and loser species indicate drivers of herb layer changes over two decades in forests of NW Germany[END_REF], very few studies have controlled for stand factors when analyzing the effect of forest continuity on biodiversity. In unfragmented forests, this may have led to premature conclusions on the relative effect of forest continuity on biodiversity.

Conclusion

In their review, Hermy and Verhenyen (2007) stated that the impact of current management activities on the distribution of ancient-forest species remained unclear. Through a dedicated sampling design, we showed that i) the effect of forest continuity was mediated by stand maturity, and ii) that local factors are extremely important in explaining biodiversity patterns. Specifically, we emphasize that, in the Northern Alps, the colonization credit of most plant and springtail species has been paid off and that, as a result, the transient biodiversity deficit usually related to forest continuity has moved toward equilibrium between ancient and recent forests [START_REF] Jackson | Balancing biodiversity in a changing environment: extinction debt, immigration credit and species turnover[END_REF]. Therefore, in unfragmented forest landscapes, above-and below-ground species richness and composition clearly appear to be governed mostly by niche processes (see also [START_REF] Jamoneau | Fragmentation alters beta-diversity patterns of habitat specialists within forest metacommunities[END_REF]. This emphasizes the necessity to better control for soil and stand factors when analyzing the response of biodiversity to forest continuity.

Finally, though it has been shown that fragmentation enhanced dispersal limitation (Honnay et al. 2002b;[START_REF] Jamoneau | Fragmentation alters beta-diversity patterns of habitat specialists within forest metacommunities[END_REF], most studies documenting the effect of forest continuity on biodiversity have used sampling designs where recent forests were disconnected from ancient forests [START_REF] Verheyen | An integrated analysis of the effects of past land use on forest herb colonization at the landscape scale[END_REF][START_REF] Kolb | Effects of environment, habitat configuration and forest continuity on the distribution of forest plant species[END_REF][START_REF] Sciama | Floristic and ecological differences between recent and ancient forests growing on non-acidic soils[END_REF][START_REF] Naaf | Colonization credit of post-agricultural forest patches in NE Germany remains 130-230years after reforestation[END_REF]. However, in cases where the extinction debt has not yet been paid off, it may be difficult to disentangle the effect of fragmentation from the effect of forest continuity per se. We show that in unfragmented forests, the legacy effect of past land-use is difficult to detect, at least for patterns of species richness and composition. Caution must be exercised when drawing conclusions about the effect of forest continuity based on the interpretation of plant assemblage shifts in similar ecological conditions, i.e. lowland fragmented forests. Despite the fact that one third of Europe's land area is mountainous and that 41 % of that area is covered by forests (EEA 2010), mountain forests have rarely been studied through the lens of historical ecology. Ecological conditions similar to the ones in our study are likely to occur in the Carpathians, the Central European Mountains, the Balkans, the Pyrenees and other areas in the Alps (e.g. [START_REF] Kozak | Forest Cover Change in the Western Carpathians in the Past 180 Years: A Case Study in the Orawa Region in Poland[END_REF][START_REF] Gellrich | Agricultural land abandonment and natural forest re-growth in the Swiss mountains: A spatially explicit economic analysis[END_REF][START_REF] Dittrich | Separating forest continuity from tree age effects on plant diversity in the ground and epiphyte vegetation of a Central European mountain spruce forest[END_REF]). More research is therefore needed to better understand the effect of forest continuity on biodiversity patterns in unfragmented mountain forests.

(

  Prop_Forest, measured within a 500 m radius) -were measured in surrounding landscape at each sampling site. Stand variables -i.e. mean canopy openness (Canop), the proportion of conifers (Tree_Compo) and the number of diameter classes (5-cm) for standing trees (Tree_Struc) -were extrapolated from stand measurements within each 20 m radius plot. Soil variables -i.e. the humus index[START_REF] Ponge | Humus Index as an indicator of forest stand and soil properties[END_REF]

FrenchFig. 2 .

 2 Fig. 2. Relative influence of parameters related to soil, stand and/or landscape on plant and springtail richness in the French pre-Alps. Sums of AICc weights are provided for 716 models based exclusively on (i) soil and/or stand variables; (ii) soil and/or landscape variables; (iii) stand and/or landscape variables. 717

Fig. 3 .

 3 Fig. 3. Constrained canonical analysis of principal coordinates of plant and springtail species composition in 718 relation to the interaction between forest continuity and stand maturity. 719

in soil, stand and landscape variables in relation to forest continuity and stand maturity

  geophytes, hemicryptophytes and chamaephytes). The sequencing of springtail metabarcodes generated 16 394 927 reads, which, after filtering, resulted in a total of 188 MOTUs (mean = 42; range 24-59) (see Appendix S1).

	Variations Soil, stand and landscape variables varied according to forest continuity and stand maturity
	classification (see Appendix S3

, 195 plant species were recorded on the 70 sites. Because overstory is directly influenced by silviculture, we focused on understory species (S = 157; mean = 27; range 8-46), i.e. herbaceous plants (therophytes,

Table 2 .

 2 The 5 among 50 top-ranking models (see Appendix S6) predicting plant and springtail richness in the

				Ancient	Recent		Mature	Overmature
	Variables		Description			p-value			p-value
				Mean (±SD)	Mean (±SD)		Mean (±SD)	Mean (±SD)
	Soil variables
	H_index	Humus index (based on humus form)	2.68 (±1.36)	2.91 (±1.68)	0.7146	3.19 (±1.39)	2.33 (±1.53)	0.0006
	pH	pH-H	2 O (topsoil)	6.31 (±1.03)	6.40 (±0.88)	0.8580	6.45 (±1.08)	6.24 (±0.81)	0.3799
	CN_ratio	Total carbon / total nitrogen ratio (topsoil)	17.49 (±2.34)	17.71 (±2.70)	0.0453	18.12 (±2.46)	17.00 (±2.45)	0.9617
	Stand variables
	Canop	Mean canopy openness (%)	10.47 (±5.24)	10.27 (±5.46)	0.7436	8.65 (±3.82)	12.31 (±6.09)	0.0061
	Tree_Compo Coniferous proportion	75.77 (±14.19) 75.51 (±17.82) 0.2164 72.37 (±16.29) 79.33 (±14.80) 0.0063
	Tree_Struc	Number of diameter classes of standing trees	11.97 (±1.76)	12.61 (±2.09)	0.5805	11.08 (±1.01)	13.61 (±1.85)	0.0001
	Landscape variables
	Alpha_LU	Land-use diversity (500 m radius)	2.51 (±1.04)	3.03 (±0.98)	0.0317	2.89 (±1.10)	2.61 (±0.97)	0.5744
	Dist_Forest Distance to forest edge in meters	226 (±111.77)	139 (±86.31)	0.0043 174 (±105.00) 197 (±113.71) 0.4129
	Prop_Forest Forest proportion (500 m radius)	92.32 (±10.81) 87.00 (±13.30) 0.0043 86.30 (±13.99) 93.74 (±8.57)	0.3190

ANOVAs, type III SS).

overmature stands in the French pre-Alps.

provided for models based exclusively on (i) soil and/or stand variables; (ii) soil and/or landscape variables; (iii) stand and/or landscape variables.

Acknowledgements

We thank Gilles Favier, Sophie Labonne, Fanny Lebagousse and Pascal Tardif from IRSTEA for help in the field. We are grateful to Delphine Rioux and Amélie Saillard from LECA for DNA extraction and to Frédéric Boyer from LECA for his assistance in bioinformation processing. We also thank Vicki Moore for correcting the English manuscript. Financial support was provided by IRSTEA, by the LTSER Zone Atelier Alpes (a member of the ILTER-Europe Network) and by Labex OSUG@2020 (ANR10 LABX56).

Data accessibility

Data can be accessed through the Dryad Data Repository: http://dx.doi.org/10.5061/dryad.6c886 [START_REF] Janssen | Data from: Present conditions may mediate the legacy effect of past land-use changes on species richness and composition of above-and below-ground assemblages[END_REF]).

Authors' contributions P.J., M.F., J.J.B. and C.B. conceived and designed the experiments. P.J., S.B., and P.T. performed the experiments. P.J. analyzed the data and wrote the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Supporting information

Additional supporting information may be found in the online version of this article: Appendix S1. Springtail assessment by environmental DNA metabarcoding. As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials may be re-organized for online delivery, but are not copy-edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Appendix

Plant richness


Eugenio Marinelli
email: eugenio.marinelli@eurecom.freddy.ghabach@eurecom.fr

Eddy Ghabach

Yiqing Yan
email: yiqing.yan@eurecom.frraja.appuswamy@eurecom.fr

Thomas Bolbroe

Omer Sella
email: osella@imperial.ac.ukt.heinis@imperial.ac.uk

Thomas Heinis

Raja Appuswamy

Digital Preservation with Synthetic DNA

Keywords: DNA storage, long-term archival, preservation, SIARD-DK

The growing adoption of AI and data analytics in various sectors has resulted in digital preservation emerging as a cross-sectoral problem that affects everyone from data-driven enterprises to memory institutions alike. As all contemporary storage media suffer from fundamental density and durability limitations, researchers have started investigating new media that can offer high-density, long-term preservation of digital data. Synthetic Deoxyribo Nucleic Acid (DNA) is one such medium that has received a lot of attention recently. In this paper, we provide an overview of the ongoing collaboration between the European Union-funded, Future and Emerging Technologies project OligoArchive and the Danish National Archive in preserving culturally important digital data with synthetic DNA. In doing so, we highlight the challenges involved using DNA for long-term preservation, and present a holistic data storage pipeline that brings together several novel techniques (standardized file storage, motif-based DNA encoding, scalable read consensus to name a few) to provide reliable, passive, obsolescence-free digital preservation using synthetic DNA.

Introduction

Today, we live in an increasingly digital society. Digital data pervades all disciplines and has established itself as the bed rock that drives our society, from enabling data-driven decisions based on machine learning, to encoding our collective knowledge compactly in a collection of bits. Thus, preservation of digital data has emerged as an important problem that must be addressed by not just memory institutions today, but also by institutions in several other sectors.

In order to preserve digital data, it is necessary to first store the data safely over a long time frame. Historically, this task has been complicated due to several issues associated with digital storage media. All current media technologies suffer from density scaling limitations resulting in storage capacity improving at a much slower rate than the rate of data growth. For instance, Hard Disk Drive (HDD) and magnetic tape capacity is improving only 16-33% annually, which is much lower than the 60% growth rate of data [START_REF] Fontana | Mooreâs law realities for recording systems and memory storage components: Hdd, tape, nand, and optical[END_REF]. All current media also suffer from media decay that can cause data loss due to silent data corruption, and have very limited lifetime compared to the requirements of digital preservation. For instance, HDD and tape have a lifetime of 5-20 years. A recent survey by the Storage and Networking Industry Association stated that several enterprises regularly archive data for much longer time frames [19]. Thus, the current solution for preserving data involves constantly migrating data every few years to deal with device failures and technology upgrades. A recent article summarized the financial impact of such media obsolescence on the movie industry [START_REF] Perlmutter | The lost picture show[END_REF].

In project OligoArchive [START_REF] Appuswamy | OligoArchive: Using DNA in the DBMS storage hierarchy[END_REF], we are exploring a radically new storage media that has received a lot of attention recently-Deoxyribo Nucleic Acid (DNA) [START_REF] Church | Next-Generation Digital Information Storage in DNA[END_REF][START_REF] Goldman | Toward Practical High-capacity Low-maintenance Storage of Digital Information in Synthesised DNA[END_REF][START_REF] Erlich | DNA Fountain enables a robust and efficient storage architecture[END_REF][START_REF] Organick | Random access in large-scale DNA data storage[END_REF]. DNA is a macro-molecule that is composed of smaller molecules called nucleotides(nt). There are four types of nucleotides: Adenine (A), Cytosine (C), Guanine (G), and Thymine (T). DNA used for data storage is typically a singlestranded sequence of these nucleotides, also referred to as an oligonucleotide (oligo). DNA possesses several key advantages over current storage media. First, it is an extremely dense three-dimensional storage medium with a capacity of storing 1 Exabyte/mm 3 which is eight orders of magnitude higher than magnetic tape, the densest medium available today [START_REF] Corporation | semiconductor synthetic biology roadmap[END_REF]. Second, DNA is very durable and can last millennia in a cold, dry, dark environment. A recent project that attempted to resurrect the Woolly Mammoth using DNA extracted from permafrost fossils that are 5000 years old is testament to the durability of DNA even under adverse conditions [START_REF] Shapiro | Mammoth 2.0: will genome engineering resurrect extinct species?[END_REF]. Thus, data once stored in DNA can be left untouched without repeated migration to deal with technology upgrades. Third, as long as there is life on earth, we will always have the necessity and ability to sequence and read genomes, be it for assembling the genome of a previouslyunknown species, or for sequencing the genome to detect diseases causing variations. As a result, unlike contemporary storage technologies, where the media that stores data and the technology to read data are tightly interlinked, DNA decouples media (biological molecules) from read technology (sequencing), thus reducing media obsolescence issues.

In this work, we provide an overview of the ongoing collaboration between the Danish National Archive and project OligoArchive in demonstrating a holistic solution for long-term preservation of culturally significant data using DNA. We present a motivating use case for long-term digital preservation, outline the challenges involved in using DNA as a digital storage medium, and present the end-to-end pipeline we have put in place to overcome these challenges.

Context and Background

Danish National Archive Use Case

The Danish National Archives is a knowledge center documenting the historical development of the Danish society. The archive collects, preserves and provide access to original data with the purpose of supporting current and future possible needs of the Danish community -public authorities as well as private citizens. A huge part of this work includes the preservation of digitally created and retrodigitized data securely and cost-effectively. Thus, the archive has received and preserved such data since the 1970s.

The archival material used for this work consists of selected hand-drawings made by the Danish king Christian IV (1577-1648). Although his reign was marked by military defeat and economic decline, Christian IV stands out as one of the most prominent, popular and admired characters in the line of Danish kings. The hand-drawings date to the period 1583-1591 where the king was 6-14 years old. The material is a part of a larger archival unit consisting of numerous documents and records4 . The specific image5 used for this experiment presents a naval battle between several warships. Besides emphasizing the young king's admiration for warfare and naval tactics, the material further indicates his high level of cultural education as well as his talent for drawing. At Danish National Archive, the material is thus ranked as having "Enestående National Betydning" (meaning unique national significance).

DNA storage challenges

Using DNA as a digital storage medium requires mapping digital data from its binary form into a sequence of nucleotides using an encoding algorithm. Once encoded, the nucleotide sequence is used to synthesize DNA using a chemical process that assembles the DNA one nucleotide at a time. Data stored in DNA is read back by sequencing the DNA molecules and decoding the information back to the original digital data.

A simple way to convert bits into nucleotides is to adopt a direct mapping that converts 2 bits into a nucleotide, for instance 00 to A, 01 to T, 10 to G, and 11 to C. This way a binary sequence is translated to an arbitrary sequence of nucleotides. However, such a simple approach is not feasible due to several biological limitations imposed by DNA synthesis and sequencing steps. First, DNA synthesis limits the size of an oligo between hundred to few thousands of nucleotides. Therefore, data must be divided into several pieces, with each piece being stored in an oligo. However, unlike current storage devices, oligos do not have logical addressing. Hence, indexing information that can help to identify the order in which the oligos, and hence the corresponding data bits, must be reassembled back during recovery must be stored together with the data bits and integrated in each oligo.

Second, oligos with repeat sequences (like ACACAC), or long consecutive repeats of the same nucleotide (like AAAAA), and oligos with extreme GC content, where the ratio of Gs and Cs in the oligo is less than 30% or more than 70%, are known to be difficult to synthesize, sequence, and process correctly. Thus, when constructing oligos, constraints need to be enforced to minimize homopolymer repeats and balance GC content. Further, care must be taken to minimize similarity across oligos as having too many oligos with positionally-similar nucleotide sequences can exacerbate sequencing errors and make it difficult to identify the original oligo.

Third, sequencing and synthesis are not error free even for well-formed oligos, as they introduce substitution errors, where a wrong nucleotide is reported, or indel errors, where spurious nucleotides are inserted or deleted. Both sequencing and synthesis also introduce bias. Some oligos are copied multiple times during synthesis, while others are not. Similarly, some oligos are read thousands of times during sequencing while others are not sequenced at all. Thus, it is important to use error correction codes in order to recover data back despite these errors.

In addition to the aforementioned media-level challenges in using DNA as a digital storage medium, there are also other problems associated with digital preservation that DNA does not solve. Any digital file stored on DNA is an encoded stream of bits whose interpretation makes sense only in the context of the application used to render, manipulate, and interact with that file format. While DNA might be able to store data for millennia, the associated applications and file formats might become obsolete. Thus, in addition to preserving data, it is also necessary to preserve the meaning of data by ensuring that data is stored in a preservation-friendly, non-proprietary format. Digital data can also be altered due to a variety of reasons and additional data-integrity techniques should be put in place to ensure that data retrieved from DNA can be trusted to be the same as the original source. The digital preservation community has long pioneered file formats, information systems, and operational methodologies for solving such format obsolescence issues [START_REF]Digital Preservation Handbook[END_REF]. Thus, a holistic DNA-based preservation solution should build on such techniques to solve both media and format obsolescence issues.

Design

In this section, we will describe the end-to-end pipeline we have put in place to overcome the aforementioned challenges.

Overcoming format obsolescence with SIARD-DK

In Denmark, all public institutions and organizations that produce data worthy of persevering are legally bound to submit them to a public archive. As the vast majority of data in the Danish public sector are organized as databases with or without files in various formats, the focus has been on archiving these data in a standardized, system-independent and cost efficient manner. As a result, the archive has implemented a Danish version of the SIARD format (Software Independent Archiving of Relational Databases) [START_REF] Of Congress | SIARD (Software Independent Archiving of Relational Databases) Version 1[END_REF] named SIARD-DK for storing of such data. SIARD is an open format, designed for archiving relational databases in a vendor-neutral form and is used in the CEF building block "eArchiving".

The first step in preserving data is extracting it and creating an SIARD-DK Archival information Package (AIP). In the creation of this particular AIP, the digitized material was converted to TIFF format. Information relevant to the images such as the preservation format of the files, their title, creator and original size, descriptive information, etc., was extracted and packaged together with relevant documentation in the AIP-format. This was done using proprietary tools developed at the Danish National Archive. The usage of SIARD for storing the files guarantee that the material is preserved in a rich format with relevant metadata stored in a standardized, system and vendor independent way. The resulting AIP is a single ZIP64 file that internally contains the TIFF images, in addition to XML and XSD files that store the schema of the archive and metadata information. This allows for strict validation of the AIP. Further, an MD5 value of each file is stored inside the archive and serves as the fixity to verify data integrity on retrieval.

DNA data storage pipeline

The end-to-end DNA media storage pipeline is presented in Figure 1. In the rest of this section, we will provide an overview of both the write path that takes as input the SIARD zip file and stores it in DNA, and the read path that restores back the zip file from DNA. Write Path In order to store the archive on synthetic DNA, the zip file is first encoded from binary into a quaternary sequence of oligonucleotides, and then synthesized to generate synthetic DNA. The steps for encoding the SIARD archive file into oligos is presented in Figure 2. During encoding, the file is read as a stream of bits and pseudo randomized. In other encoding methods, randomization is used as a way to limit the number of homopolymer repeats in

Consensus Sequencing

Oligoarchive

Zip file

Stream of bits

Randomized bits LDPC encoding Indexing Oligos Constrained code Fig. 2. Encoding bits into oligos each oligo. In our encoding, homopolymer repeats are handled by an inner constrained code that we explain later. Thus, we do not need randomization for avoiding homopolymer repeats. We use randomization primarily to improve the accuracy of our clustering and consensus methods in the data decoding stage.

As mentioned before, data stored in DNA is read back by sequencing the DNA to produce reads, which are noisy copies of the original oligos that can contain insertion, deletion, or substitution errors. Our read clustering and consensus methods rely on the fact that the original oligos are well separated in terms of edit distance so that the distance between a noisy read and its corresponding oligo is much smaller than the distance between two oligos. This assumption makes it possible to cluster similar reads and infer the original oligo with a high accuracy. A long sequence of zero or one bits can violate this assumption as they can lead to multiple oligos being similar, or even identical, to each other. A long sequence of bits can also lead to oligos with repetitive sequences (example: ACA-CACAC). This can pose problems during data decoding, especially if paired-end sequencing is used, where the DNA is partially read from either direction (5" to 3" and 3" to 5") with overlap. In such a case, the two reads corresponding to each orientation must be merged into a single representative read. Repetitive sequences can create issues during this merging process. While we can develop more advanced solutions to perform merging, randomization provides a simple solution that eliminates such issues while ensuring that similarity across oligos is also minimized. After randomization, error correction encoding is applied to protect the data against errors. We use large-block length Low-Density Parity Check (LDPC) codes [START_REF] Gallager | Low-density parity-check codes[END_REF] with a block size of 256,000 bits as the error correction code, as it has been shown to be able to recover data in the presence of intra-oligo errors, or even if entire oligos are missing [START_REF] Chandak | Improved read/write cost tradeoff in dnabased data storage using ldpc codes[END_REF]. We configure LDPC to add 10% redundancy to convert each sequence of 256,000 bits into 281,600 bits with data and parity. Each 281,600 bit sequence is then used to generate a set of 300-bit sequences, where each 300-bits is composed of 281 data bits and a 19-bit index that is used to order the sequences. Each 300-bit sequence is then passed to a constrained code that converts it into an oligonucleotide sequence.

The constrained code essentially views each oligo as a concatenation of several shorter oligonucleotide sequences, that we henceforth refer to as motifs. In our current configuration, the constrained code breaks up each 300-bit sequence into a series of ten 30-bit integers. Each 30-bit integer is fed as input to a motif generator that takes the 30-bit value and produces a valid 16nt (nucleotides) motif as output. The motif generator does this mapping by pre-constructing an associative array where the 30-bit value is the key and a 16nt motif is the value. This array is built by first enumerating all possible motifs of length 16nt. Then, all motifs that fail to meet a given set of biological constraints are eliminated.

Our current motif generator is configured to allow up to two homopolymer repeats (AA,CC,GG, or TT), and admits motifs with a G-C content in the range 0.25 to 0.75. With these constraints, using 16nt motifs, out of 4 16 possible motifs, we end up with only 1,405,798,178 unique, valid motifs with which we can encode any possible 30-bits of data ( as the number of all possible 30-bits values is 2 30 = 1, 073, 741, 824 which is lower than 1,405,798,178. Thus, we use the billion motifs as values in the array corresponding to keys in the range 0 to 2 30 . The 30-bit value is thus used as the key into this associative array to produce the corresponding 16nt motif. Thus, at the motif level, the encoding density is 1.875 bits/nt. The reason we limited ourselves to 16nt and 30bits is the fact that this associative array occupies around 100GB of memory, which we can easily meet using our current hardware. The motif generator can be extended to larger motif sizes and more relaxed biological constraints which can lead to higher bit densities. But as this would require the use of external storage, we leave this open to future work.

Using the associative array, each 300-bit sequence is encoded as concatenation of ten motifs, each with a length of 16nt, leading to an oligo that is 160 nucleotides long. We would like to explicitly point out here that the length of an oligo is a configurable parameter. Thus, while we use 160nts in our current system due to favorable pricing provided by our synthesis provider, our encoder can generate shorter or longer oligos if necessary, and automatically adjust various aspects (like the 19-bit index and 281-bit data size) based on desired oligo length.

Read Path To retrieve back the SIARD archive, the DNA is sequenced in order to retrieve back the nucleotide sequence of oligos. As mentioned before, sequencing produces reads, which are noisy copies of the original oligos. Thus, we need a consensus procedure to infer the original oligos from the reads. In prior work, we structured this process as sequence of three algorithms, as depicted in Figure 3. First, we identify all pairs of strings that are similar to each other. As modern sequencers produce hundreds of millions of reads, this first task is extremely computationally intensive due to use of the edit distance as a metric for comparing strings. Thus, we have developed an efficient similarity join algorithm, called OneJoin [START_REF] Marinelli | Onejoin: Cross-architecture, scalable edit similarity join for dna data storage using oneapi[END_REF], that exploits the fact that due to randomization during encoding, reads corresponding to the same original oligo are "close" to each other despite some errors and "far" from the reads related to other oligos. The results obtained from the join algorithm are then used to quickly identify clusters of strings that are similar to each other. Each cluster thus groups all reads belonging the same oligo. Finally, we apply a position-wise consensus procedure that uses multiple reads to infer the original oligo in each cluster. While our prior approach was able to infer original oligos with a high accuracy, there were two problems. First, we found that under some datasets, particularly for high coverage reads, OneJoin's memory and computational usage were too high. As OneJoin is a string similarity join, it produces as output all possible pairs of reads that are similar to each other. Thus, given a coverage N , the computational and memory requirements of OneJoin were O(N 2 ). Second, the use of a general purpose string similarity join led to functionality repetition at multiple places in the read path. For instance, OneJoin internally uses an edit distance check to filter out strings that are not similar. Later in the read pipeline, we had to repeat the edit distance computation in the alignment stage once to get position-wise consensus. This repetition led to needless overhead. To solve these problems, we have developed a new consensus procedure that we refer to as OneConsensus.

In the following sections, we describe the key stages of OneConsensus algorithm. In order to efficiently identify and group all the similar reads, our algorithm relies on two well-known algorithmic tools that allow to drastically cut down the computational time: CGK-Embedding and Locality Sensitive Hashing (LSH).

CGK Embedding. As we mentioned in the previous section, the similarity metric used in OneConsensus is the edit distance. Given two strings x and y, the edit distance is defined as the minimum number of edit operations i.e. insertions, deletions and substitutions, necessary to transform x in y. Another metric, commonly used to compare strings is the Hamming distance. However, the latter takes into account only the number of mismatches between the two strings, or in other words the number of substitutions to transform x in y. For example, given the two strings ACACT and GACAC, their Hamming distance is 5 since there are no matches, but the edit distance is 2 since it suffices to add G and remove T

From these definitions, we can make the observation that the edit distance takes into account information about the ordering of characters and captures the best alignment between two strings. However, while Hamming distance has complexity that is linear with the string length, edit distance has complexity that is quadratic. While several dynamic programming optimizations exist for accelerating edit distance computations [START_REF] Ukkonen | Algorithms for approximate string matching[END_REF], they often rely on pre-specified distance thresholds and are unable to provide performance competitive with fine-tuned Hamming distance computations in practice. Given the complexity of this metric, we rely on randomized embedding techniques to minimize the overhead of the edit distance computations.

Randomized embedding refers to a set of methods that map a complex metric space into a simpler one. CGK-embedding algorithm, recently proposed by Chakraborthy et al. [START_REF] Chakraborty | Streaming algorithms for embedding and computing edit distance in the low distance regime[END_REF] is one such algorithm that can map problems from an edit space into a Hamming space. Given two strings x, y of length N taken from an alphabet such that d E (x, y), the edit distance between x and y, is less than K, CGK-embedding is a function f :

N → 3N that maps strings x and y into f (x) and f (y) such that, with probability at least 0.99, the Hamming distance of d H (f (x), f (y)) is bounded by K 2 when d E (x, y) < K. This implies that the distortion D, defined as the ratio D(x, y) = d H (f (x),f (y))

d E (x,y)
, is at most K. Thus, as long as the edit distance is small, the distortion of embedding is small. This implies that the Hamming distance of embedded strings will accurately track the edit distance of the original strings, thereby making it possible to replace the expensive edit distance computation with cheap Hamming distance computation [START_REF] Zhang | Embedjoin: Efficient edit similarity joins via embeddings[END_REF]. end if 8:

Algorithm 1 CGK-embedding

i ← i + Rj 9: end for 10: return S ′ The pseudo-code of embedding algorithm is shown in Algorithm 1. In this case the procedure is applied to all strings of length N that are composed of the characters A, C, G, T , representing the DNA alphabet. Given an input string, the algorithm builds the corresponding embedded representation by appending one character at time taken from the input string. The character appended can be the repetition of the previous character or the next character in the input string according to the value of a binary random string. In other words, the pointer of the current character in the input string increases or remains the same depending on the random string value, that can be 0 or 1. When the pointer to the input string goes beyond the string length, the embedded string is padded with a special character P . In general, P can be any character that is not included in the alphabet of the given dataset S. For sake of simplicity, in Algorithm 1 we use 0 for padding. In essence, what we get as the output of embedding is a string where characters from the input string can be repeated one or more times.

LSH for Hamming Distance. One of the main advantages of moving from the edit distance space to the Hamming space is being able to use some useful algorithms that are valid in the Hamming space only and instead not applicable to the edit distance. One of these algorithms is Locality Sensitive Hashing (LSH) [START_REF] Gionis | Similarity search in high dimensions via hashing[END_REF]. Definition 1. We call a family H of functions (d 1 , d 2 , p 1 , p 2 )-sensitive for a distance function D if for any p, q ∈ U (where U is the item universe):

-if D(p, q) ≤ d 1 then P [h(p) = h(q)] ≥ p 1 ,
that is, if p and q are close, the probability of a hash collision is high;

-if D(p, q) ≥ d 2 then P [h(p) = h(q)] ≤ p 2 ,
that is, if p and q are far, the probability of a hash collision is low;

where h ∈ r H are hash functions randomly sampled from the family of hash functions H

Considering two bit-string p and q of length N. In the Hamming distance case, the hash function is defined as the i th bit of these strings. Thus, if their Hamming distance is d H (p, q), that is, the number of bits that differ position wise in the two strings, then the probability that any given bit at a random position is the same in both strings is 1 -d H (p,q) N . Thus, the bit-sampling LSH family for Hamming distance, defined as:

H N = {h i : h i (b 1 ...b N )) = b i | i ∈ [N ]} is d 1 , d 2 , 1 -d1
N , 1 -d2 N -sensitive for the two Hamming distances d 1 < d 2 . We use Hamming LSH over embedded reads to separate out the reads into different buckets such that with a very high probability, reads within a bucket are similar to each other, and hence correspond to the same reference.

Clustering based on Edit Distance At this stage, we have all reads grouped in hash buckets based on their similarity. However, we still have two problems to solve. (1) LSH can produce false positives, meaning that two dissimilar reads can end up in the same bucket. The main consequences is that if reads are very different, the consensus procedure lead to the wrong result. ( 2) Reads can have different lengths due to insertions and deletions errors. Thus, we need to adjust the reads in order to make their lengths uniform while taking into account possible insertion/deletion errors. Both these problems can be solved by aligning the reads in each bucket. More specifically, given a bucket, we sort the reads based on length such that reads with length matching the reference oligos are moved to the front of the bucket. Then, starting with the first read in the bucket, we align all the following reads to the first one. The intuition behind sorting the reads is as follows. If a read has the same length as the reference oligos (160nt), either the read has no insertion/deletion errors, or there are an even number of insertion and deletion errors. Since the probability of errors is low with shortread sequencing, the former scenario is more likely. Thus, by picking reads that are of the correct length and aligning the rest of the reads to it, we increase the probability of finding the correct original reference oligo.

The alignment of reads gives two pieces of information: the edit distance and the Compact Idiosyncratic Gapped Alignment Report (CIGAR). The edit distance allows us to identify reads that are actually dissimilar even if they are in the same buckets. We group only similar reads to form a cluster. The CIGAR contains the base-by-base alignment information (the sequence of matches, insertions and deletions) needed to align one read to the other. Using this information, we adjust the reads by adding gaps where there is a deletion error, or deleting nucleotides where there is an insertion error. Once all similar reads are found within the bucket, we save the cluster and remove the reads from the bucket. Any dissimilar reads that were not a part of the cluster are still left in the bucket. In order to deal with false positives by LSH, the remaining reads are then processed again in the same way, with the procedure being repeated until the bucket is empty.

Position-Wise Consensus

The result of the previous stages is a set of clusters. All reads within a cluster are noisy copies of the same reference oligo with some errors in random positions. At this point of the algorithm the only missing step is the consensus procedure. The consensus algorithm works on a per position basis. The key idea is that despite some random errors, all reads in a cluster are aligned and adjusted based on the edit distance and CIGAR. Thus, if we consider a specific position across all reads, it is likely that the majority of reads in that position will contain the correct nucleotide while only some of them report the wrong one. This implies that for any given position, it will be enough to take the most-frequent nucleotide as the consensus outcome. Repeating this procedure for each position, we produce one inferred oligo per cluster. We would like to point out here that not all oligos need to be correctly inferred. In fact, as we show later, some original oligos might not appear at all in the inferred set, and other inferred oligos might have errors. We rely on the parity added by LDPC codes at a higher level to recover data despite these errors.

The inferred oligos are then passed to the decoder which reverses the steps shown in Figure 2. The constrained code is first used to convert each 160nt oligo back into 300-bit sequences by converting each individual 16-nt motif into its corresponding 30-bit value. The index stored in each 300 bits is used to reassemble bits back in the correct order. The LDPC decoder is then used to recover back data even if some bits were wrongly decoded, or some bits were zeroed out as corresponding oligos were missing. The decoded data is then derandomized to obtain a stream of bits that corresponds to the original input.

Evaluation

At this stage of our collaboration, we have assembled the entire pipeline. We are in the process of carrying out a real-life, large-scale synthesis experiment. As we do not have results from real experiments yet, we will provide a preliminary, simulation-driven evaluation in this section. Note that we simulate only the synthesis and sequencing steps in Figure 1. We encode/decode the real dataset using our pipeline.

Experimental Setup

The raw SIARD archive that is fed as input to our pipeline is 12.9MB in size. With redundancy added by LDPC, the resulting binary data to be stored on DNA is 15.19MB in size. We encode the SIARD archive generating 405,212 oligos, each with a length of 160nts. Using these original oligos, we then generate four million reads by using a short-read simulator6 tool, that adds random errors such as insertion, substitution, and deletion in each read to mimic the actions of an Illumina DNA sequencer. This corresponds to an average coverage of 10×, meaning that each oligo, on average is covered by 10 noisy copies.

Benchmark with sequence alignment

We begin our analysis by visualizing the coverage across oligos. While the average coverage is 10×, the overall coverage typically follows a negative binomial distribution, with some oligos being covered hundreds of times, and some being not covered at all. To visualize this, we aligned the reads to reference with BWA-MEM v0.7.17 [START_REF] Li | Aligning sequence reads, clone sequences and assembly contigs with bwamem[END_REF], a state-of-the-art short-read aligner. Based on the alignment result from BWA-MEM, we show the histogram of coverage across oligos in Figure 4. The x-axis is the coverage (or number of reads that map to an oligo), and y-axis is the number of oligos with that coverage in log scale. As can be seen, the coverage distribution spans a range from 1 to 26, with majority of oligos being covered 5-15× as expected given that the simulator was configured to produce 10× coverage.

We also use the alignment result to show the histogram of number of errors (indels and substitutions) in Figure 5, where x-axis represents the number of errors, and y-axis represents the number of oligos with that error count in log scale. As can be seen, the error distribution is right skewed, with a majority of reads having fewer than 2 mismatches, and a few reads having as many as 8 mismatches which accounts for a 5% error rate given the oligo length of 160. Finally, we also used the alignment result to verify that each read uniquely maps to an oligo (absence of 'XA' field in the SAM output file produced by BWA-MEM). This shows that oligos are "far" from each other in terms of edit distance due to randomization. BWA-MEM performs base-by-base, edit-distance-based base-by-base alignment and reports CIGAR and alignment As a result, it is computationally very intensive. While such an alignment is required for genomic data analysis in order to determine the exact location of a read in the genome, for the purpose of DNA storage, we found that it is sufficient to simply map each read to an oligo without full alignment. Accel-Align (v1.1.1; [START_REF] Yan | Accel-align: A fast sequence mapper and aligner based on the seed-embed-extend method[END_REF]) is a short-read aligner that we have developed in the context of project OligoArchive that supports alignment-free mapping-only mode that can quickly map reads to oligos. Thus, we present a comparative analysis of BWA-MEM and Accel-Align here with the goal of presenting it as a open-source tool that can be used by other researchers for both DNA storage and more broadly, for analyzing genomic data. We use the BWA-MEM and Accel-Align to align/map reads to oligos. Using BWA-MEM's alignment as the gold standard, we evaluate Accel-Align's accuracy. Table 1 shows the execution time and the percent of correctly aligned reads. It shows that Accel-Align can performing mapping 6x faster than BWA-MEM at a slight drop in accuracy of 0.004%. On futher analysis, we found this drop to be due to Accel-Align's inability to map reads with very high error rate. Accel-Align is specifically built to trade off performance and accuracy for Illumina-based short-read sequencing reads whose error rate is typically much lower than 5%. As our simulated reads offer a more pessimistic error model, Accel-Align experiences a slight drop in accuracy. Despite this, as the coverage histogram with Accel-Align and BWA-MEM are near identical, we have found Accel-Align to be a very useful tool for analyzing both DNA storage reads and more broadly, genomic data [START_REF] Yan | Accel-align: A fast sequence mapper and aligner based on the seed-embed-extend method[END_REF].

End-to-end Decoding Results

Having presented an analysis of the reads, we will now present the decoding results. Using OneConsensus procedure described earlier, we obtain the inferred oligos using the simulated dataset. Table 2 shows the error statistics for these oligos. The figures are obtained by comparing the inferred oligo for a certain index with the corresponding original reference oligo. We see that we are able to infer 404,075 oligos that correspond to 99.7% of the original oligos perfectly without errors. In addition, 1010 oligos were inferred with some errors and 127 oligos were completely missing. Note that errors in an oligo does not imply that the entire oligo is different from the original, but differs only with respect to a few motifs. For this reason, we also report the difference between inferred data and original encoded file in terms of number of bits.

We then use the constrained code to convert these inferred oligos into 300-bit sequences, and reassemble them in order based on the 19-bit index. At this stage, we will certainly have situations where an oligo is missing due to sequencing simulation bias, or an oligo could not be converted back into 300-bits due to errors. In both cases, there will be a corresponding index whose data bits cannot be recovered. We insert a sequence of zero bits for such indices and use the reconstructed binary together with the LDPC decoder to restore the original input. Despite the errors reported in Table 2, the LDPC decoder was able to recover back the original archive completely, thanks to the additional parity information added during encoding. Finally, in this section we compare OneConsensus with OneJoin-based consensus algorithm. Table 3 compares the two algorithms in terms of the accuracy achieved in inferring the encoding oligos. In terms of the number of oligos correctly inferred in their entirety, meaning an exact match between the inferred oligo and the original reference oligo, we see that OneConsensus slightly underforms OneJoin. But as we mentioned earlier, oligos can also differ by just a few motifs only, making the statistic about the correctly inferred oligos insufficient to determine the actual accuracy of the two algorithms. For this reason, we report also the number of missing oligos and the number of bits wrongly inferred by the two algorithms. The number shows that OneConsensus outperforms OneJoin, as it mistakes only misses 127 oligos, while OneJoin 595 oligos. Overall, OneConsensus leads to 42678 bit errors, compared to the 95935 bit errors produced by the OneJoin-based consensus.

In terms of memory consumption, we observed that OneJoin reaches a peak of 2.5GB, while OneConsensus requires only 1.1GB. While the difference may not seem too striking for this dataset, we would like to point out that while OneConsensus memory consumption grows linearly with the dataset size, One-Join requires an amount of memory that is quadratic with coverage.

Conclusion and Future Work

In this work, we provided an overview of the ongoing collaboration between project OligoArchive and the Danish National Archive in using DNA to preserve culturally significant digital data. Building on prior work on molecular information storage and digital preservation, we presented a holistic, end-to-end pipeline for preserving both data and the meaning of data on DNA, and tested the pipeline using simulation studies. There are several avenues of future work we are pursuing. First, as described earlier, we are in the process of carrying out a large-scale experiment to validate our pipeline using real data. Second, we are investigating various optimizations to both encoding and consensus algorithms to support alternate synthesis and sequencing technologies with potentially higher error rates. Finally, while we addressed the question of preserving data in this work, we left open the question of preserving the decoding algorithm itself. Recent work has investigated the design of nested universal emulators that can be used to preserve and emulate such decoders using analog media like film or archival paper [START_REF] Appuswamy | Universal layout emulation for long-term database archival[END_REF]. Thus, we are investigating methods to combine DNA-based digital data storage with analog media-based decoding logic storage.

Fig. 1 .

 1 Fig. 1. DNA Storage Pipeline

Fig. 3 .

 3 Fig. 3. Various steps in the OneJoin consensus procedure.

3N 1 : i ← 0 2 :

 12 Input: A string S ∈ {A, C, G, T } N , a random string R ∈ {0, 1} 3N and a char for padding P = 0 Output: The embedded string S ′ ∈ for j = 0 → 3N -

Fig. 4 .

 4 Fig. 4. Histogram of occurrence of inferred oligos

Fig. 5 .

 5 Fig. 5. Histogram of mismatches

Table 1 .

 1 Performance and accuracy of BWA-MEM and Accel-Align(map mode)

		Exec. time (second) Correctly aligned (%)
	BWA-MEM	25.5	100
	Accel-Align	3.8	99.9

Table 2 .

 2 Statistics for decoding of SIARD archive

	#Original Oligos	405212
	#Correctly Inferred Oligos 404075
	#Incorrectly Inferred Oligos 1010
	#Missing Oligos	127
	#Incorrect bits	42678

Table 3 .

 3 Statistics for OneConsensus and OneJoin-based consensus

		OneConsensus OneJoin
	#Correctly Inferred Oligos	404075 404104
	#Incorrectly Inferred Oligos	1010	513
	#Missing Oligos	127	595
	#Incorrect bits	42678 95935

https://www.flickr.com/photos/statensarkiver

https://www.flickr.com/photos/statensarkiver/28273082238

https://sourceforge.net/projects/bbmap/

Acknowledgments

This work was partially funded by the European Union's Horizon 2020 research and innovation programme, project OligoArchive (grant agreement No 863320).