dataesr/scientific-paragraphs-categorization huggingface dataset:
As rarefaction can introduce sequencing data bias (McMurdie and Holmes, 2014), we addressed this potential issue using a combined strategy to assess the effect of sequencing depth on sampling completeness and on the dispersion induced on diversity estimation. First, we performed bootstraps on rarefaction curves to observe ASV counts by increasing the library size to see if a plateau was reached (Fig. A3 ). Then, using the Mirlyn package  (Cameron et al. , 2021)  we performed stepwise bootstraps on read numbers per sample to assess Shannon diversity index value changes by increasing sequencing depth (Fig. A4 ). Finally, in each individual sample we tested the dispersion induced by the rarefaction process at different sequencing depth (Fig. A4 ). All of these elements were used to determine the rarefaction level.
Network construction Before network construction, we subsetted new ASV tables from the original dataset for each sampling site (Nat 1, Nat 2, OF, LR, EB) and compartment (rhizosphere and root) to avoid confusing environmental filtering effects. We used rhizosphere samples for network construction as the soil-plant root interface is known to be the most responsive to environmental changes compared to roots  (Shi et al. , 2016)  we also noted this feature in our dataset. Rare taxa were eliminated (<10 reads across samples) and prevalence thresholds of 20 % for 16S and 15 % for ITS were used. Networks were constructed using SPIEC-EASI  (Kurtz et al. , 2015) , which is robust against compositional data by using clrtransformation to remove the unit-sum constraint of compositional data, as proposed by Aitchison  (Aitchison, 1981).
Acknowledgement We would like to thank all those involved in the SoFunLand project, and the RCU ( Royal Commission for AlUla)  and Afalula ( Agence Française pour le Developpement d'AlUla)  for their financial support. We would also like to thank  O. Domergue ,  R. Destabel  and  A. Dehail  for their help with soil and biomolecular analyses. And our invaluable local guide  Abdulaziz Marzouq Salman AlNajem  for his knowledge of the field. Finally, we would like to thank  A. Camuel  for her advice on the manuscript and  D. Manley  for his English proofing.
Data availability Demultiplexed sequence data are available at NCBI (Bio-project PRJNA1028866, Biosamples SAMN37854590 to SAMN37854744). Soil XRF raw datas are available in Supplementary file 1. Extended PER-MANOVA are available in Supplementary file 2. Phylum level of positive and negative interactions proportion for each bacteria network are listed in Supplementary file 3. Phylum level of positive and negative interactions proportion for each bacteria-fungi networks are listed in Supplementary file 4. ANCOVA results for the stability analysis of each network with the paired comparisons with the estimated marginal means method is available in Supplementary file 5. Module and network hubs taxonomical information, topology, Z i and P i values are listed in Supplementary file 6. Proportion of taxa roles in each network is available in Supplementary file 7.
Image registration The evaluation of bimodal tomography data sets generally requires image registration. This involves adjusting the spatial resolution and precisely aligning the data to obtain congruent three-dimensional X-ray and neutron images of the samples. Registration of bimodal data sets is particularly challenging because the gray scale and contrast levels of individual phases can vary substantially from one modality to another. We used the registration software called "spam"  (Stamati et al. , 2020)  to register the set of samples containing MP particles in sand enriched with organic materials. The remaining samples, i. e. , the soil samples containing MP film fragments, MP fibers and lupine roots, were registered using the software "Fiji". The field of view, spatial resolution, and 3D orientation of the volume datasets were manually adjusted similar to the procedure described in  Haber-Pohlmeier et al.
(2019). The X-ray tomography (CT) data was used as reference during image registration with "spam", which means that the voxel size, the canvas size, and the number of tomographic slices of the neutron tomogram were adapted to the dimensions of the CT images. The multimodal registration algorithm is considered as an extension of classical digital image correlation methods  (Roubin et al. , 2019). Its goal is to determine the linear transformation operator F, affecting the coordinates of one modality image, in our case neutron tomography (NT), so that it matches the reference image, here X-ray CT. Prior to the actual registration in "spam" the NT data were scaled to the spatial resolution of the CT dataset and roughly aligned manually using "Fiji". For robustness and speed, "spam" prefers to register highly downscaled images, which eliminates high-frequency variations  (Roubin et al.
MP segmentation and 3D visualization Data visualization and segmentation of MP particles were performed using the rendering software "VGSTUDIO MAX" (Volume Graphics, Heidelberg, Germany, version 3. 1). In terms of 3D image processing, the program allows versatile user interaction directly at the level of the reconstructed 3D sample volume. A 3D rendered view of the reconstructed sample volume is generated by setting various rendering parameters, such as gray-scale value thresholds and opacity curves in the 3D histogram. This proved to be very beneficial when using the implemented region-growing algorithm for segmentation, which is initiated by manually positioning the growth seeds. To enable reliable subjective identification of MPs, a 3D view was generated showing only the components with the strongest neutron attenuation (c.
Detection of microplastic particles in root-soil systems Finally, we tested the detection of differently shaped microplastic particles, specifically polycarbonate fragments and PP film fragments, in root soil systems. For this purpose, ten-day-old lupine plants grown in a mixture of sandy soil and MP particles were imaged at the low soil moisture level (0. 04 cm 3 cm -3 ), which allowed good contrast between MP and soil and between roots and soil. A particular challenge in segmenting this neutron image data is distinguishing the root structures from the MP, since the root tissue and the plastic material have similar attenuation properties (Fig. 7a , for the polycarbonate MP). First, the roots were segmented in the neutron image (Fig. 7b ). The root system has a clearly defined tap root and a large number of lateral roots, some of which have started to grow along the container walls.
In this context, various sources of Earth observation data can provide information with great potential to characterise coastal morphological changes over large areas and with high temporal resolution. This is the case of the optical imagery of the Sentinel-2A and 2B satellites (supplied by the European Commission through the European Space Agency, ESA), as well as the Landsat series (supplied by NASA). All of them acquire images covering large extensions (i. e. , tens of kilometres) with short revisit times (5 days for Sentinel-2 when combined at the equator, with shorter revisit intervals at higher latitudes  (Bergsma and Almar, 2020). Large SDS packages may allow to characterise the morphological response to storm events and anthropogenic actions  (Cabezas-Rabadán et al. , 2019a; b; Pardo-Pascual et al.
, 2014)  and their consequences on the beach functions  (Cabezas-Rabadán et al. , 2019c). The Shoreline Analysis and Extraction Tool  (SAET, Palomar-Vázquez et al. , 2023)  has been developed to carry out the automatic detection of the shoreline positions associated with the storm events using the images from the Sentinel satellites of the EU Copernicus program and the Copernicus Contributing Missions. Based on the tool SHOREX  (Cabezas-Rabadán et al. , 2021; Sánchez-García et al. , 2020)  SAET has been developed within the framework of the ECFAS project (A proof of concept for the implementation of a European Copernicus coastal flood awareness system, GA n • 101,004,211, www. ecfas. eu) focusing on achieving high autonomy and efficiency, capacity to be applied on large coastal segments, and robustness when working on different coastal types.
, 2006). The aforementioned limitations cannot be avoided, and every accuracy evaluation must deal with them. The aim of the present work is to evaluate the automatically extracted SDSs using the new software SAET. The results are focused on (a) the effect of the input images TOA and BOA and (b) the influence of the initial pixel-level shoreline on the final shoreline. The latter effect is influenced by both the band or index, the thresholding method, and the morphological filter used for deriving the pixel-level shoreline from the land/water mask. The tests are run over different sites that enable the analysis of the robustness of the resulting shorelines over different environments. The work is organized in different sections. Firstly, the sites in which the assessment was carried out are presented together with the very high resolution (VHR) images employed as reference data.
Definition of SDSs The satellite-derived shorelines were defined from the freely available optical Sentinel-2 imagery using the software SAET. This tool provided the framework for downloading the source S2 images as well as the algorithm to obtain the SDSs. The workflow followed four main phases (Fig. 3  For the SDSs assessment, the SDS extraction was carried out by combining different S2 processing levels and extraction methods (Fig. 4 ). Just like the reference shorelines, the SDSs were defined as the instantaneous water/land interface. The workflow in SAET followed 4 steps.
This observation agrees with the evaluation performed by  Bishop-Taylor, et al. (2019)  and the satisfactory performance in the extraction of large shoreline datasets in intermittent shallow lakes using SHOREX  (Palomar-Vázquez et al. , 2022). Nonetheless, AWEInsh was expected to offer good results as it was designed to highlight the radiometric differences between the sea and the land  (Feyisa et al. , 2014)  by giving more weight to the value of the SWIR band. This way, AWEInsh would provide an improvement over the MNDWI index previously used with satisfactory results by other tools such as CoastSat  (Vos et al. , 2019a). AWEInsh with a constant value of 0 as proposed by  Feyisa et al. (2014)  makes it possible to perform the automatic thresholding required for defining the location of the APS.
The CLMS coastal zone database is still under validation, therefore some beaches are not properly mapped or a few  beach polygons are misplaced In the case of the Mira site, a better definition of the boundaries of the beach resulted in the removal of the wrong points and to the variation of the error from -133 ± 93 m to 28. 7 ± 14. 6 m (Fig. 9 ). The results are in line with  Hagenaars et al. (2018)  as the shoreline extraction is challenging in beaches with a very gentle slope (as the sites 1-3). In The Hague, the worst SDS accuracy appears when the shorelines are obtained using the AWEIsh index over the TOA images by dilating the water mask. These results are unexpected as AWEIsh is used in the presence of shadows that are not present in the analysed images. Additional validation is needed to understand if intertidal morphologies, saturated sand and puddles are treated by the algorithm as shadows.
Acknowledgements This research has been supported by the funds of the following projects:  ECFAS  (A proof-of-concept for the  implementation of a European Copernicus Coastal Flood Awareness System)  project (htt ps://www. ecfas. eu/), which has received funding from the  EU H2020 research and innovation programme  under grant agreement no. 101004211 ;  MONOBESAT  ( PID 2019-111435RB-I00 ) funded by the  Spanish Ministry of Science, Innovation and Universities ; the  Margarita Salas  contracts within the  Re-qualification programme  by the  Spanish Ministry of Universities  also financed by the  EU -NextGenerationEU , and the grant  Primeros Proyectos de Investigación  ( PAID-06-22 ),  Vicerrectorado de Investigación de la Universitat Politècnica de València (UPV). Authors acknowledge  Mercator Ocean International  and  Puertos del Estado  for providing tidal and wave conditions, and  ESA  for the satellite imagery within the framework of ECFAS.
[4]. Technical replicates (n = 3) were mixed in equimolar amounts for sequencing. Libraries were sequenced using Illumina MiSeq platform with paired-end v3 chemistry (300 bp paired-end reads). Read-pair overlapping and initial sequence data filtering was carried out using vsearch  [59]  and clustering into sequence variants (SVs) using SWARM v2  [44] , followed by removal of likely chimeric sequences and singletons, and finally curation using LULU  [25] , as previously described by  Lanzén et al. (2020). Curated SVs were classified using CREST  [40]  v3. 2. 2 with SilvaModPR2 v138 as reference database (https://github. com/lanzen/CREST). The resulting SV contingency tables were further processed in R using the package vegan v2. 5  [53]. Relative abundances were calculated using the vegan function decostand (with mode="total").
Identification of Hg-pollution indicators To identify potential indicator taxa of Hg pollution in the 16 S rRNA gene amplicon sequencing dataset, we used Threshold Indicator Taxa Analysis as implemented in TITAN2 v2. 1  [5]. Average relative abundances were used for this analysis to pool biological replicates lacking replicated measurements of Hg concentration. This analysis was based on a custom Mercury Pressure Index (Hg PI) derived from measurements of total Hg concentrations of each sample and transformed according to the ecotoxicological Threshold Effect Limit (TEL) of Hg in sediments (0. 13 ppm as established by the National Oceanic and Atmospheric Administration  [50]. For details on pressure index calculation see  [8]. Taxa significantly associated with high Hg PI (p < 0. 01) in all samples were considered as potential Hg-pollution bioindicator taxa.
Screening of merA, merB and hgcAB genes in metagenomic contigs and taxonomic affiliation Six sediment samples collected in winter 2013 from area of study (EN17, EN20, EOK10, LOK10, EOI20 and EU08) and characterized by different levels of environmental impact  [39] , were used for metagenomic analysis. Metagenomes were obtained by shotgun sequencing with 150 bp paired-end sequencing on the Illumina HiSeq Platform, in the National Center for Genomic Analysis (CNAG) in Barcelona. Metagenomic reads were quality controlled using FastQC  [1] , and sequencing adaptors trimmed using Cutadapt v1. 15  [45]. Resulting sequence data was subjected to serial co-assembly as described in Tully et al. [68]. Briefly, sample-wise assemblies were first carried out using Megahit  [42]. All resulting contigs across assemblies with a size of 2000 bp or larger were then concatenated and subjected to co-assembly using Minimus2  [63].
To determine the normalized relative abundance of each contig (Reads Per Kilobase per Million mapped reads, RPKM), all reads were mapped to contigs using Bowtie2 v. 2. 1. 0  [38]  and resulting counts divided by the total number of mapped reads from each dataset, as well as the length of the corresponding contig. Open reading frames (ORFs) were predicted from all resulting coassembled contigs (including "singleton" contigs, i. e. original samplespecific contigs) and translated to amino acid sequences using Prodigal  [32]. The resulting ORFs were screened for merA, merB and hgcAB genes. Custom Hidden Markov Models (HMM) for merA and merB were built querying UniProtKB for the KEGG orthologs (KO) identifier K00520 and K00221, respectively (BioStudies accession number S-BSST1154). Confirmed sequences for each gene (377 and 55, respectively) were aligned in mafft v7.
402 with default parameters  [37]. In the case of hgc genes, we used the procedure described in Capo et al. [13]. Briefly, HMM profiles of hgcA and hgcB genes derived from the Hg-MATE database v1. 01142021  [27]  were used to develop a HMM using hmmbuild from HMMER v3. 2. 1  [24]. We considered genes with E-values < 10 -3 as significant hits. To decrease false positives, we then used the stringency cutoff defined by Capo et al. [13] , to verify the presence of the conserved motifs from hgcA (NV(I)WCA(A/G/S)GK) and hgcB genes (C(M/I)EC(G/S)(A/G)C), and thereby performed a manual inspection of the presence of hgcA and hgcB genes in our dataset. Certain hgcA genes were found side-by-side with hgcB genes on the same contig. Each putative merA and merB gene retrieved from the metagenomes was taxonomically assigned by mapping them to UniRef90 (release 2021_03 from 9 of June 2021), with mmseqs2 development version, commit 13-45111, with the taxonomy workflow options "-max-accept 100 -tax-lineage 1 -e 1E-5 -v 3 -a" and converted to table with mmseqs createtsv.
All ranks out of domain, phylum, class, order, family, genus or species were stripped out from classification and missing fields were marked as "unclassified". The lowest common ancestor for each sequence was also recorded. In the case of hgc genes, phylogenetic analyses were done with the pplacer approach as part of the marky-coco pipeline  [13]. The software MEGA11  [66]  was used for peforming the alignment of merA genes extracted from the metagenomes and amplified from the isolates (see above) using the UPGMA method, and for building the phylogenetic tree (ML method).
Association networks of Desulfobacterota taxa For reconstruction of association networks of Desulfobacterota in the sediment samples, we used FlashWeave v1. 4  [64]  with default parameters on the 16S rRNA gene amplicon sequencing dataset. The network was constructed based on OTUs resulting from merging all SVs with the same taxonomic grouping. The redox potential and total Hg concentration were included as parameters and submitted to FlashWeave in addition to OTU abundance tables. Associations with p>0. 01% were discarded. Negative associations (edges) were also discarded to simplify the network structure and minimize potential bias introduced by differences in community composition associated with contrasting environmental conditions unrelated to Hg concentration. The resulting networks were inspected and illustrated using Cytoscape v3.
However, hgcB was not found side-by-side to hgcA in this MAG. The other two MAGs (EN20maxbin. 003 and LOK10maxbin. 002) were also assigned to Desulfobacterota (Desulfuromonadia and Synthrophobacteria, respectively). One of these MAGs, MAG EN20maxbin. 03, contained merA and hgcAB genes, indicating a potential to carry out both mercury reduction and methylation (Supplementary Table  S3 ). This MAG featured 10. 05% contamination as estimated by the software checkM, based on the Fig. 3. Normalized abundances (RPKM) of merA, merB and hgcA-carrying contigs in each sample (A) and its taxonomic assignment (B), retrieved from metagenomic analyses in the six samples analyzed. The gradient of mercury concentration in the samples (represented by a triangle in pannel A) ranged from 108 to 2100 ng/g in dry weight (see Supplementary Table  S1 ).
Acknowledgements This work has been co-funded by the  European Climate, Infrastructure and Environment Executive Agency  ( CINEA ,  EMFF-BlueEconomy-2018-863584 , project  MER -CLUB). We received additional support from  Basque Government  (projects microgAMBI and  BIOMATRIX23 ), the  Basque Monitoring network  financed through a collaboration between the  Basque Water Agency (URA)  and  AZTI  and  EU-MINECRITICAL  ( C17. 17 ,  MET2021-00-000 ). A. L. is supported by  IKERBASQUE (Basque Foundation for Science). We thank  Bastien Duval  for his support in Hg analysis and  Andrea G. Bravo  for useful discussions on the results of the investigation. This is contribution number 1194 of AZTI,  Marine Research, Basque Research and   Technology Alliance (BRTA).
Data Availability We have shared the data/code in the manuscript with their accession numbers. Sequence reads and derived MAGs are available from the Sequence Read Archive of the INSDBC via EBI WebIn (BioProject accessions PRJEB65319 and PRJEB62120). HMM profiles for merA and merB genes can be found in BioStudies under accession number S-BSST1154. 3. merA sequences obtained by Sanger sequencing have been deposited in NCBI under accession numbers OR672480-OR672517.
The northern permafrost region has been projected to shift from a net sink to a net source of carbon under global warming. However, estimates of the contemporary net greenhouse gas (GHG) balance and budgets of the permafrost region remain highly uncertain. Here we construct the first comprehensive bottom-up budgets of CO2, CH4, and N2O across the terrestrial permafrost region using databases of more than 1000 in-situ flux measurements and a land cover-based ecosystem flux upscaling approach for the period 2000-2020. Estimates indicate that the permafrost region emitted a mean annual flux of 0. 36 (-620, 652) Tg CO2-C y-1, 38 (21, 53) Tg CH4-C y-1, and 0. 62 (0. 03, 1. 2) Tg N2O-N y-1 to the atmosphere throughout the period. While the region was a net source of CH4 and N2O, the CO2 budget was near neutral with large uncertainties.
Study area The spatial extent of permafrost defined in this study includes areas within the northern permafrost region as defined in  Obu et al. (2021)  and restricted to the Boreal Arctic Wetlands and Lakes Dataset area that had the key land cover classes for our upscaling  (BAWLD, Olefeldt et al. 2021a,b)  (Fig. 1 ). As a consequence, the BAWLD-RECCAP2 permafrost region does not include large areas underlain by permafrost in Central Asia and the Tibetan plateau. The BAWLD-RECCAP2 permafrost region considered in this study is 18. 5 million km 2. All flux estimates were scaled to the BAWLD-RECCAP2 permafrost region (hereafter permafrost region). The study area overlaps several other RECCAP2 regions  (Ciais et al. 2022) , and no specific effort to harmonise the budgets presented here with the RECCAP2 budgets of those regions is made in this paper.
GHG budgets from ecosystem flux upscaling Data-driven ecosystem flux upscaling of GHG budgets for a reference time period of 2000-2020 was calculated by summing up flux budgets from terrestrial ecosystems, inland waters, lateral fluxes, fire emissions, and geological fluxes. To calculate the total net regional GHG flux (F x ), we used the following equation: 𝐹 = 𝐴 𝑗 × 𝐹 𝑗𝑥 where F x is annual permafrost region gas flux for the GHG species of interest x, A j is the area of each land cover class j (Fig. 2 , Table  S1 ), and F jx is the land cover average GHG flux density for species x (Table  S1 ). We used existing synthesis databases and upscaled gridded data products published in the past five years to estimate annual and growing season mean fluxes per land cover type. All budget numbers are presented as the weight of C and N (i.
GHG fluxes from terrestrial land cover types The land cover classification used for the analysis was adapted from the Boreal-Arctic Wetland and Lake Dataset (BAWLD) land cover  (Olefeldt et al. 2021a,b). The BAWLD land cover classes are distinguished based on moisture regime, nutrient/pH regime, organic-soil depth, hydrodynamics, and the presence or absence of permafrost  (Olefeldt et al. 2021a ). To match the observational GHG flux datasets, we simplified the nine terrestrial land cover classes in BAWLD into five: Boreal Forests; Non-permafrost Wetlands; Dry Tundra; Tundra Wetlands; and Permafrost Bogs (Fig. 2 ). Classes were defined as: • Non-permafrost Wetlands include permafrost free bogs, fens, and marshes with no nearsurface permafrost (see Canadian Wetland Classification system). (Schiller and Hastie, 1996; Simpson et al.