dataesr/scientific-paragraphs-categorization huggingface dataset:
2016). Although the emergence of predominant scales is significant in our dataset, it does not explain a large proportion of the variation found (Fig. 3 ). In this regard, it is important to note that many regular patterns found were only apparent when windowing the subplots studied, which highlights the importance of selecting relevant sizes of the observation windows in different situations (e. g. depending on the plant form that is aimed to be studied). Moreover, this suggests that regular and irregular patterns are often mixed in the field (something already suggested by  Tarnita et al. 2017). This finding calls for future studies on the fact that multi-specific sites in real ecosystems may show both regular and irregular spatial patterns, which is seldom considered but may have important consequences to determine the validity of spatial patterns as tools to monitor ecosystem degradation.
Indeed, some studies have linked patch formation with increases in the ability of plants to maintain a high biomass due to an increased resource capture efficiency  (Aguiar and Sala 1999 , Boer and Puigdefábregas 2005 , Meyra et al. 2012 ) and regular spotted patterns have been interpreted in some simulations as alternative states of desert systems  (Rietkerk et al. 2004 , Kéfi et al. 2010). Although aridity per se was not a significant predictor of the shape of PSDs, it strongly modulated the effects of biotic attributes on the shape of PSDs in shrublands, but not in grasslands (Fig. 4 ). Grasslands in our dataset were particularly prevalent under moderate aridity conditions (Supplementary material Appendix 4 Fig. A8 ), so the lack of interactions with aridity might be a consequence of a smaller aridity range in the grasslands surveyed.
Whereas grasses might depend strongly on water pulses, shrubs have deeper roots and do not depend so strongly on rainfall seasonality  (Sala et al. 1997 ). On the other hand, sand content showed a positive effect on the emergence of regular spatial patterns in shrublands, and a negative effect on grasslands. Some studies report associations between loamy soils and the occurrence of regular spatial patterns which concur with those reported in grasslands  (Deblauwe et al. 2008). We think two processes may explain these contrasting results for shrublands. First, sand content may increase experienced water stress and constraint soil amelioration. Shrublands in our database comprised larger aridity gradients than grasslands, so plants may be more prone to compete than to facilitate in the most extreme environments.
Funding -This work was funded by the  European Research Council  ( ERC Grant  agreements  242658  [BIOCOM] and 647038 [BIODESERT]). MB was supported by a  FPU fellowship  from the  Spanish Ministry of Education, Culture and Sports  (ref. AP2010-0759 ). FTM acknowledges support from a  Research Award  from the  Alexander Von Humboldt Foundation. SS was supported by the  Spanish Government  under a  Ramón y Cajal contract  ( RYC-2016-20604 ). The research of SK has received funding from the  European Union 's  Seventh Framework Programme  ( FP7/2007-2013)  under grant agreement no. 283068  (CASCADE).
Topography data and analysis The topography data were taken from the global digital elevation model GTopo30 (https://lta. cr. usgs. gov/GTOPO30, last access: 20 August 2018). GTopo30 has a horizontal grid spacing of 30 arcsec corresponding to a spatial resolution between 928 m in the East/West direction at the equator, 598 m at WHI (50 • N) and 373 m at the SUM polar station (72. 6 • N). In the North/South direction, 30 arcsec are almost constant with latitude and correspond to 921 m at the equator and 931 m at the poles. The geographical coordinate system WGS84 (World Geodetic System revised in 1984) from GTopo30 was projected in the Universal Transverse Mercator (UTM) conformal projection to ensure homogeneity in vertical and horizontal coordinates. Due to the altitude averaging over each grid cell, there is typically an altitude difference between the true station altitude and its corresponding grid location.
The mean and median of the differences between the station altitude and its representative grid point are 190 m (8. 6 %) and 140 m (5. 8 %). For stations situated near the summit, the difference can be significant: five stations have an underestimation of their altitude larger than 500 m corresponding to 15 %-32 % (see the Supplement and particularly Table  S1  for further details), despite, according to the manual, the high reported GTOPO30 accuracy (minimum accuracy of 250 m at 90 % confidence level with a RMSE of 152 m). The TopoToolbox-master version of the free shareware TopoToolbox (https://topotoolbox. wordpress. com/, last ac-cess: 20 August 2018), which is a set of matlab functions offering analytical GIS utilities in a non-GIS environment  (Schwanghart and Scherler, 2014) , was used as a principal tool for the topographic relief and flow pathways analysis in the Digital Elevation Model (DEM) analysis.
Aerosol parameters Aerosol datasets from 25 high altitude and 3 mid-altitude stations (Table  1 ) were available for this study, 21 of them coming from GAW (Global Atmospheric Watch) stations. The datasets comprise absorption coefficient, scattering coefficient and/or number concentration and cover time periods ranging from at least a year up to more than a decade of measurement (see Table  S3  in the Supplement). Stations with time series shorter than one year were not used, as their data are not representative of a complete seasonal cycle. Due to the non-normal distribution of the aerosol parameters, the 5th, 50th and 95th percentiles were taken as representative of the minimum, central position and maximum concentrations. No correction for standard pressure and temperature was applied in order to use the measured aerosol properties and concentration at high altitude.
For consistency, the measured hourly absorption and scattering coefficients were adjusted to a wavelength of 550 nm if reported at a different wavelength using an assumed scattering or absorption Ångström exponent of 1. Additionally, the scattering coefficients were corrected for angular truncation error. Due to the measurement technique and the low aerosol concentrations at many high altitudes, filter-based photometers regularly measure negative absorption coefficients at some of these sites. Some datasets contain up to 20 %-30 % of negative absorption values. Depending on the data owner's policy, these negative values were either left in the dataset, set to zero (or to a minimum value) or considered as missing values. To ensure a similar treatment for all datasets, negatives values, zeros and minimum values attributed to negatives were thus set to missing values.
The diurnal and seasonal cycles were only analyzed on datasets longer than 2 years. To be able to statistically calculate the diurnal and seasonal cycles, the autocorrelations at 1 h (first lag) were first removed from the dataset by a whiten-ing procedure  (Wang and Swail, 2001). The autocorrelations at each lag time were then calculated on the whitened dataset taking into account missing data (see Supplement for further explanations). Only autocorrelation values statistically significant at 95 % confidence level were kept. As the diurnal (24 h) and annual cycles (365 days) were not well defined due to variable meteorological conditions and some shorter datasets, the autocorrelation at lags of 22 to 26 h and at lags of 350 to 380 days were summed to obtain the strength (i. e. , the cycle amplitude) of the diurnal and seasonal cycles, respectively.
As NPF can be enhanced at low temperature, (i. e. , NPF seasonal cycle can be anticorrelated with the seasonality of the CBL height) the relationship between number concentration (and, to a lesser extent, scattering) with diurnal cycle strength may be obscured by seasonal changes in atmospheric dynamics. The only SS correlation with station altitude is found for the scattering coefficient seasonal cycle. Similar to the anticorrelation of G8 with the aerosol parameter percentiles, there is also a high anticorrelation between the particle number concentration diurnal cycles and G8 suggesting that the local slope steepness inhibits both the transport of polluted air masses and NPF. Apart from a correlation at 90 % confidence level between DBinv and the absorption coefficient, the lack of further SS correlations for topographic parameters with the aerosol parameter seasonal cycles can be attributed to several factors: (i) the relatively small time period (2-5 years) covered by most of the datasets leading to difficulties in the statistical determination of a yearly periodicity due to inter-annual variability, (ii) the low aerosol concentration at high altitude sites often resulting in measurements near the detection limits of the instruments (see for example the problem with the absorption coefficient at Sect.
nilu. no (last access: August 2018). The OMP and SBO data will be soon accessible on the WDCA web page, but data can be also obtained by contacting Paolo Fiahlo (fialho. paulo@gmail. com) and Gehrard Schauer (gerhard. schauer@zamg. ac. at), respectively. PEV data will be soon be accessible online on the Bolin Centre database (https://bolin. su. se/data/, last access, August 2018) and can be obtained by contacting Radovan Krejci (radovan. krejci@aces. su. se). MUK data can be obtained by contacting Hooda Rakesh (rakesh. k. hooda@fmi. fi). Supplement. The supplement related to this article is available online at: https://doi. org/10. 5194/acp-18-12289-2018-supplement. Author contributions. DA, MA, HA, NB, ME, PF, HF, AGH, RH, IK, RK, NHL, AM, JM, NAN, MP, VP, LR, SR, GS, KS, SS, JS, PT and PV contributed to the measurement of the aerosol time series.
Swedish participation was supported by the Swedish Foundation for International Cooperation in Research and Higher Education (STINT) and the Swedish research council FORMAS. -CMN: the European Commission funded ACTRIS, ACTRIS-2 EU project and NEXTDATA National project funded by MIUR. -IZO: Global Atmospheric Watch program funded by AEMET and by the project AEROATLAN (CGL2015-66299-P), funded by the Ministry of Economy and Competitiveness of Spain and the European Regional Development Fund (ERDF). - -MUK: Ministry of Foreign Affairs of Finland, Academy of Finland (project no. 264242) for the financial support and an esteemed collaboration of FMI and TERI. -NCOS: Max Planck Institute for Chemistry and the Chinese Academy of Sciences (Project SKLCS-ZZ-2017 and KJZD-EW-G03-03). -PDI: Federal Office of Meteorology and Climatology Me-teoSwiss through the project Capacity Building and Twinning for Climate Observing Systems (CATCOS) Phase 2, Contract no.
JFJ: International Foundation High Altitude Research Stations Jungfraujoch and Gornergrat (HFSJG) and the Federal Office of Meteorology and Climatology, MeteoSwiss, within the Swiss program of the Global Atmosphere Watch (GAW) of the World Meteorological Organization, funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 654109 (ACTRIS2), and the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract no. 15. 0159-1. The opinions expressed and arguments employed herein do not necessarily reflect the official views of the Swiss Government. -LLN: the Taiwan Environmental Protection Administration. -MSA and MSY: MINECO (Spanish Ministry of Economy and Competitiveness), the MAGRAMA (Spanish Ministry of Agriculture, Food and Environment), the Generalitat de Catalunya (AGAUR 2014 SGR33 and the DGQA) and FEDER funds under the PRISMA project (CGL2012-39623-C02/00).
This work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 654109. Marco Pandolfi is funded by a Ramón y Cajal Fellowship (RYC-2013-14036) awarded by the Spanish Ministry of Economy and Competitiveness. -MLO, NWR and SUM: the NOAA observatory staff at each station and the NOAA Climate Program Office's Atmospheric Chemistry, Carbon Cycle and Climate (AC4) program.
81025332  between the  Swiss Agency for Development and Cooperation and MeteoSwiss ,  Günther Wehrle  from  Paul Scherrer  institute, the site operator  Luis Valle  from  DGAC  and to the maintenance staff from MeteoChile. -WLG:  National Scientific Foundation of China  (  41675129 ),  National Key Project of MOST  ( 2014CB441201 ) and the  CMA Innovation Team for Haze-fog Observation and Forecast. -ZEP: Support from  Swedish Environmental Protection Agency (Naturvårdsverket) ,  The Swedish research council FORMAS, Knut and Alice Wallenberg Foundation (KWA)  and  Norwegian Polar Institute (NPI) -ZUG:  H. E. Scheel  for performing aerosol measurements at ZUG and  Wolfgang Junkermann ,  Karlsruhe Institute of Technology, Garmisch-Partenkirchen, Germany , who give us access to this dataset. Edited by: Laurens Ganzeveld Reviewed by: two anonymous referees.
De novo metatranscriptome assembly and differential expression Trimmomatic  [36]  was used to trim nucleotides of poor quality and reads < 40 bp were removed. Reads from all 12 biological samples were assembled into a de novo transcriptome using Trinity software set to default parameters  [3]. Assembled contigs shorter than 200 bp were discarded. Bowtie2  [37, 38]  was used to align the RNA-seq reads back to the de novo transcriptome with -a -X 600 parameters on top of default parameters. Low count contigs, below 29 total counts across all libraries, were filtered out (corresponding to those contigs below the 80th percentile  [39] ). Read back-mapping rates were an average of 64. 91% across all 12 root samples. Raw and normalised transcript abundance (tpm) was calculated using eXpress  [40]  with default parameters.
Differential expression (DE) was estimated using EBSeq based on median normalisation  [41]  and on an empirical Bayes model framework  [42, 43] , keeping contigs with a posterior probability of differential expression (PPDE) ≥ 0. 95 (with target false discovery rate controlled at 5%). EdgeR  [44]  was used to generate MA plots based on trimmed mean of M-values normalisation method  [45]. Extended quality control, assembly and normalisation information is provided in Additional file 1, and scripts are provided in Additional file 2.
Annotation The metatranscriptomic (unconstrained) annotation strategy, which queries a broad range of protein sequence repositories, was performed as outlined by Gonzalez et al. [46]. Briefly, the de novo assembled contigs were annotated using three major protein databases (nr NCBI non-redundant protein database, SwissProt and TrEMBL) as well as the S. purpurea 94006 reference genome. During the informed annotation step, Populus trichocarpa reference genome was also added. UniProt Archive (UniParc) database was used to protein blast differentially expressed (DE) contigs that did not have a hit in any of the databases. NCBI nucleotide database was used to nucleotide blast all DE contigs. E values < 10 -4 (protein blast) and < 10 -6 (nucleotide blast) were used as cutoffs. A previously reported method for selecting annotation (based on the percentage of maximum potential bitscore) from blastx returns was used to help select the primary annotation given multiple high scoring alignments for a single sequence, all statistically characterised as nonrandom  [2].
Results and discussion Assembly, mapping and annotation A total of 189,849 confident contigs were confidently assembled de novo from 456,182,049 sequence reads. Overall, 125,151 contigs (66%) were confidently annotated as potential protein encoding transcripts, but 64,698 contigs (34%) had no protein homologue in SwissProt, Trembl or nr databases and are referred to as unknown henceforth (a substantial number of these were DE and investigated in greater detail in Additional files 1 and 3). The annotated transcripts included 91,053 Eukaryota (48%), 33,222 Bacteria (18%), 187 Archaea (0. 01%) and 139 viral sequences (0. 01%; including the common Illumina spike) (Fig. 1 ). Transcripts best annotated from eukaryotes were, as could be expected, mostly derived from Viridiplantae  (46, 817) , although a close second was Fungi  (40, 352) , with Metazoa (1417), Amoebozoa (918), Stramenopiles (617), Euglenozoa (390), Alveolata  (198)  and others (such as metagenomic or environmental samples) comprising only minor proportions of the transcript diversity (324).
Full annotation is provided in Additional file 4 and an interactive Krona of total annotation is available at: https://github. com/gonzalezem/Tripartite_Metatranscriptomics_article identification investigations in Additional file 1). Others have expressed even stronger scepticism, such as Bent et al. [51] , who discussed the difficulties in providing reliable diversity indices from microbial fingerprinting methodology  [52] [53] [54] [55] [56] [57]. An alternative, top-down approach is employed here which accepts uncertainty of the microbiome community, assembling as many contigs as possible within the system and performing differential expression analysis on all retained contigs without any necessity for a priori sequence information before independent annotation using the world's major protein repositories.
Transcripts clearly identified as deriving from Salix (the mass majority from the Salix purpurea 94006 genome) made up 85. 6% (40,098) of the contigs from plant species. While this was expected, it is perhaps more informative to recognise that 14. 4% of plant transcripts were best annotated outside of Salix, mostly from the very close relative Populus (9. 0%) but also from other genera (Vitaceae 0. 6%; Cenchrinae 0. 4%; Funariaceae 0. 4; Acalypheae 0. 4%; Oryza 0. 2%; Amygdaleae 0. 2%; Arabidopsis 0. 2% and others 2. 8%). This is roughly in line with previous research suggesting between 6 and 10% of Salix purpurea genes can be better identified by querying the broader protein repositories  [2, 46]  and that, importantly, over 10% of likely Salix gene expression can be obscured by only mapping to a reference genome, even when species specific.
Metaorganism patterning in gene expression A total of 12,576 contigs were identified as DE between non-contaminated and contaminated soil conditions (Fig. 2 ). DE contigs (putative transcripts when annotated as protein encoding) may be responding directly to hydrocarbon toxicity or to environmental alterations due to the high soil concentrations of petroleum hydrocarbons, such as nutrient or water availability. Beyond this, differential expression could also be due to secondary responses representing the shift in the microbiome community or, indeed, expression of other organisms independent of the shift in community within the metaorganism. No sequences were discarded based on annotation, especially since previous studies have demonstrated the potential for the removal of sequences (unexpected a priori) to biologically  [60]  and technically  [2]  confound data integrity and interpretation.
Pooling of highly complex gene expression data from non-model plants into ontology groups can give the impression of an overview of the biological system (GO terms are provided for query Additional file 4);   [47] ; the proportion of each taxonomic grouping is defined by the number of distinct contigs, whereas the colour represents the relative abundance (transcripts per million, tpm) of transcripts in each taxon. An interactive Krona graph to assist navigation of DE contig annotation origin is available at: https://github. com/gonzalezem/Tripartite_Metatranscriptomics_article. A full contig list including expression information, annotation (1°and 2°) and gene ontology is provided in Additional file 4 whereas DE only contigs are provided in Additional file 9 however, an attempt to explore DE function at a transcript level across all the organisms within the samples is presented here.
Top: fold change (FC log 2 ) distribution of DE genes contaminated (black) and non-contaminated (gold). Bottom: mean transcript counts (tpm) difference between conditions against fold change per DE contig. The highly abundant transcripts discussed within the text are labelled. A full DE transcript list including expression data, functional description (if available), gene ontology terms (if available) and secondary annotation (if available) is provided in Additional file 5. PIP plasma membrane intrinsic protein (aquaporin), GST glutathione S-transferase Glycine max  [62] , can be limiting for vacuolar loading of the isoflavone conjugates, and the corresponding conjugate-hydrolysing β-glucosidase is localised to the plant root apoplast, intriguing in terms of the extensive interaction with microbes evidently ongoing within the roots.
4b  and Additional file 6). These species both dominated fungal differential expression and had high enough levels of primary annotated contigs (2303 contigs and 1077 contigs) to merit very confident identification of these species, or of close relatives, within the system. Transcript levels were almost all more abundant within roots of noncontaminated trees (94. 1 and 98. 4% for P. omphalodes and T. melanosporum, respectively). While an alteration in gene expression in response to hydrocarbon conditions is possible, the extremity of polar expression might also suggest a lower abundance for these organisms due to the impact of contamination stress and/or increased competition from more contamination-tolerant life forms. Tuber has been previously identified as associated with Salix  [83, 84] ; nonconstitutively expressed ECM high-abundance markers were DE in Tuber and were indeed in higher abundance in non-contaminated conditions, such as RAS protein (c594475_g3_i1)  [85] , implying successful Tuber ECM interaction was underway in non-contaminated conditions but was suppressed by hydrocarbon contamination.
[86] ). An interactive chord diagram and Krona graph to assist more comprehensive navigation of taxonomy and fungal secondary annotation are available at: https://github. com/gonzalezem/Tripartite_Metatranscriptomics_article. A full fungal DE contig list including expression information, annotation (1°and 2°) and gene ontology is provided in Additional file 6 whereas a full list of Basidiomycota DE contigs upregulated in roots of contaminated trees is provided in Additional file 7 Upregulated Basidiomycota function Given that the intricacy, or perhaps futility, in differentiating between saprophytes and ECM fungi has been well discussed  [89, 90] , determining the fungal mode of action from expression study is problematic  [91]  and is further confounded by the generally accepted belief of multiple ECM evolutionary events in Basidiomycota (  [92] ).
The scale of mycorrhizospheric complexity is very high; once DE genes have been sub-selected as best annotated from fungi (keeping in mind that 2575 DE contigs were unknown, having no known sequence similarity), further sub-selected based on Basidiomycota annotation and even further sub-selected based on response to contamination (just those contigs in higher abundance in contaminated conditions), 1745 DE-annotated transcripts remain to describe potential functionality (driven by positive expression) within the system. Of these 1745 DE transcripts, 70% (1227 contigs) had unknown function and were annotated as hypothetical, predicted, putative or uncharacterised proteins from Basidiomycota species. These poor functional description terms are selected against during annotation if an equally good hit is available in secondary annotation, in practice re-mining any confident functional descriptions available within the major protein data repositories.