Benjamin Alric
email: benjamin.alric@irstea.fr

Cajo J F Ter Braak

Yves Desdevises

Hugo Lebredonchel

Stéphane Dray

Investigating microbial associations from sequencing survey data with co-correspondence analysis

Keywords: co-correspondence analysis, co-occurrence network, next-generation sequencing, microbial eukaryotes, Mamiellophyceae, Prasinovirus

Introduction

Microbial communities are highly diverse [START_REF] Rappé | The uncultured microbial majority[END_REF] and drive the major ecosystem functions (e.g., carbon sequestration, climate regulation, gas regulation, nutrient cycling; [START_REF] Ducklow | Microbial services: Challenges for microbial ecologists in a changing world[END_REF][START_REF] Falkowski | The microbial engines that drive earth's biogeochemical cycles[END_REF][START_REF] Hutchins | Microorganisms and ocean global change[END_REF]. Understanding how these systems are structured and identifying the underlying processes is a crucial task to predict communities and ecosystem responses to global change [START_REF] Fuhrman | Microbial community structure and its functional implications[END_REF][START_REF] Graham | Microbes as engines of ecosystem function: When does community structure enhance predictions of ecosystem processes[END_REF]. Biotic interactions across taxonomic groups (e.g., predation, parasitism, mutualism or competition) are of broad interest because they are expected to influence the structure and composition of communities [START_REF] Wisz | The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling[END_REF]. Unfortunately, our understanding of the underlying assemblage rules of microbial communities is still limited [START_REF] Little | Rules of engagement: interspecies interactions that regulate microbial communities[END_REF][START_REF] Cordero | Microbial interactions and community assembly at microscales[END_REF].

The emergence of high-throughput sequencing techniques (next-generation sequencing; NGS) gave access to the diversity of whole microbial communities, including the non-cultivable fraction [START_REF] Handelsman | Metagenomics: Application of genomics to uncultured microorganisms[END_REF]Zimmerman, Izard, Klatt, Zhou, & Aronson, 2014;[START_REF] Zinger | Two decades of describing the unseen majority of aquatic microbial diversity[END_REF]. With the large amount of data generated in a single NGS experiment, powerful statistical methods are needed to assess and explain structural patterns in such complex data sets [START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF][START_REF] Paliy | Application of multivariate statistical techniques in microbial ecology[END_REF]. A common approach is to combine NGS techniques with network analysis to represent and characterize interactions between partners in microbial communities [START_REF] Cardona | Network-based metabolic analysis and microbial community modeling[END_REF][START_REF] Vacher | Learning ecological network from next-generation sequencing data[END_REF]. Various computational methods have been developed to infer networks from NGS data sets (e.g., CoNet: [START_REF] Faust | Microbial co-occurrence relationships in the human microbiome[END_REF]SparCC: Friedman & Alm, 2012;REBACCA: Ban, An, & Jiang, 2015;CCLasso: Fang, Huang, Zhao, & Deng, 2015; SPIEC-EASI: [START_REF] Kurtz | Sparse and compositionally robust inference of microbial ecological networks[END_REF].

Another popular approach is to use ordination methods to extract information from NGS data and describe the variations in community composition among samples [START_REF] Paliy | Application of multivariate statistical techniques in microbial ecology[END_REF]. Ordination methods arrange objects in a multidimensional space using directly the original raw data table (e.g., principal component analysis, correspondence analysis), or after computing a distance matrix (e.g., non-metric multidimensional scaling, principal coordinate analysis) [START_REF] Legendre | Numerical ecology[END_REF].

These different approaches are based on a single table composed of read counts for each Operational Taxonomic Unit (OTU) measured for several samples.

Technological advances now make it possible to acquire NGS data on different taxonomic groups simultaneously for the same samples [START_REF] Fierer | Metagenomic and Small-Subunit rRNA analyses reveal the genetic diversity of bacteria, archaea, fungi, and viruses in soil[END_REF] and lead to analyze a pair of tables (i.e., OTUs composition for the same sampling sites for two different taxonomic groups). To analyze such pair of tables, a common practice consists in merging the two tables into a single one and then applying network analysis (e.g., [START_REF] Kueneman | Inhibitory bacteria reduce fungi on early life stages of endangered Colorado boreal toads (Anaxyrus boreas)[END_REF][START_REF] Banerjee | Netowrk analysis reveals functional redundancy and keystone taxa amongst[END_REF][START_REF] Ma | Geographic patterns of co-occurrence network topological features fro soil microbiota at continental scale in eastern China[END_REF] and/or multivariate analysis [START_REF] Cannon | Dynamic microbial populations along the Cuyahoga River[END_REF][START_REF] Bergelson | Characterizing both bacteria and fungi improves understanding of the Arabidopsis root microbiome[END_REF]. However, this data aggregation is unsuitable especially when NGS data sets, which are a function of sequencing depth [START_REF] Ni | How much metagenomic sequencing is enough to achieve a given goal[END_REF], are standardized by dividing read counts by the total number of reads in each sample. In this case, the normaliszation step and further analysis are very sensitive to the difference in number of OTUs and associated counts in each taxonomic group. Hence, it is important to use techniques allowing for the analysis of a pair of NGS data tables while preserving the original structure of the data. In addition, these approaches must be able to mitigate the statistical bias stemming from high-dimensionality (i.e., a number of samples substantially lower than the number of variables), sparsity (i.e., a high proportion of zero counts), and the compositional nature (i.e., a non-independence of relative abundances induced by the row-sum normalisation) that characterize NGS data [START_REF] Li | Microbiome, metagenomics, and high-dimensional compositional data analysis[END_REF]. For network analysis, the SPIEC-EASI method has been adapted to infer associations among microorganisms in a cross-domain analysis [START_REF] Tipton | Fungi stabilize connectivity in the lung and skin microbial ecosystems[END_REF]. For multivariate analyses, several two-table methods exist and have been presented to microbial ecologists in methodological reviews [START_REF] Ramette | Multivariate analyses in microbial ecology[END_REF][START_REF] Buttigieg | A guide to statistical analysis in microbial ecology: a community-focused, living review of multivariate data analyses[END_REF][START_REF] Paliy | Application of multivariate statistical techniques in microbial ecology[END_REF]. However, these works focused on asymmetric methods (e.g., canonical correspondence analysis and redundancy analysis) that aim to explain the composition of microbial communities by a limited number of environmental predictors. Unfortunately, these methods are not adapted to link two NGS data tables as they require that there are fewer predictor variables than samples [START_REF] Dray | Co-inertia analysis and the linking of ecological data tables[END_REF] and thus are not able to deal with the highdimensionality of NGS data. Moreover, these methods compute linear combinations of the predictor variables, which is not suitable if the table of predictors contains community data that display unimodal structure and/or are sum normalised (as the NGS data).

This study aims to propose an analytical framework based on co-correspondence analysis (CoCA;[START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF], a two-table coupling method developed in community ecology, to study the distributions and assemblages between two microbial communities. This framework is based on correspondence analysis, a method that effectively handles proportional data that contain many zeroes [START_REF] Gauch | A comparative study of reciprocal averaging and other ordination techniques[END_REF][START_REF] Jackson | Compositional data in community ecology: the paradigm or peril of proportions[END_REF][START_REF] Greenacre | Power transformation in correspondence analysis[END_REF], like NGS data [START_REF] Paulson | Differential abundance analysis for microbial marker-gene surveys[END_REF]. We show how this method allows to extract information about the co-structure among two microbial communities to estimate the congruence between them. Finally, we show that the outputs of the method can be used to reorder co-occurrence networks inferred by network analysis to enhance the visualization of microbial association and the understanding of assemblage patterns within networks. Hence, our approach echoes to the current lively debate about the practices to create network visualizations which are both aesthetically appealing and have high information content (see Pocock et al., 2016 for review). We illustrate our approach using two real data sets, one on autotrophic and heterotrophic microbial eukaryotes in shallow freshwater systems and another on microalgae and viruses in marine systems.

Materials and methods

Studying cross-taxon congruence by co-correspondence analysis

Co-correspondence analysis [START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF]) is part of the class of canonical analyses with the feature to be designed to analyze a pair of tables containing abundance data and to study the co-variations between two types of communities (e.g., plants and pollinators). CoCA is based on correspondence analysis [START_REF] Benzécri | Statistical analysis as a tool to make patterns emerge from data[END_REF][START_REF] Hill | Reciprocal averaging: an eigenvector method ordination[END_REF] and preserves its fundamental properties of weighted averaging and the use of the chi-square (χ 2 ) distance for both tables (Supporting Information). The χ 2 distance is particularly adapted to NGS data as it handles properly zero values (and in particular double absences) and thus is not hampered by zero inflation [START_REF] Legendre | Numerical ecology[END_REF].

Here, co-correspondence analysis in its predictive form (pCoCA) is used for the microalgae-virus data set while the symmetric form (sCoCA) is used for the microbial eukaryote data set. Let X =[ x ij ] and Y =[ y ik ] be n × p and n × q tables containing the relative abundances of each p and q species of two communities measured at the same n samples. In the first case, pCoCA is chosen because we wish to investigate the composition of virus communities under the hypothesis that the occurrence of viruses depends mainly on whether microalgal hosts are present or absent at a particular sample location. In the second case, sCoCA is chosen as we simply wished to study the relationships between autotrophic and heterotrophic microbial eukaryotes. In its predictive form, CoCA is based on partial least-squares regression analysis (PLS) in a SIMPLS version [START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF] to deal with high-dimensionality and subsequent collinearity in the table of explanatory variables. PLS searches for a linear regression model from a set of orthogonal components (called latent factors) built from collinear explanatory variables with the constraint that these components maximize the covariance with the response variables [START_REF] Martens | Reliable and relevant modelling of real world data: A personal account of the development of PLS Regression[END_REF]. In its symmetric form, CoCA fits in the framework of co-inertia analysis (COIA, [START_REF] Dolédec | Co-inertia analysis: an alternative method for studying species-environment relationships[END_REF] that is not affected by the problem of collinearity [START_REF] Dray | Co-inertia analysis and the linking of ecological data tables[END_REF]. Co-inertia analysis relates two data tables in a symmetric way, by providing a decomposition of the co-inertia criterion on a set of orthogonal axes on which sample scores are projected [START_REF] Dray | Co-inertia analysis and the linking of ecological data tables[END_REF].

In practice, CoCA identifies associations (or common ecological gradients) between two types of biological assemblages from the same samples, by seeking the factorial axes that maximize the covariance between the weighted average sample scores (projection of rows of one table onto the factorial axes) of one community with those of the other community (Supporting Information, Eq. 12). This requires to determine species scores (projection of columns of one table onto the factorial axes) of one table as weighted averages of sample scores of the other table and sample scores as weighted averages of the species scores of their own table (Supporting Information, Eq. 13-16). So, pCoCA and sCoCA are a weighted version of PLS and COIA, respectively.

Given that CoCA is related to the correspondence analysis, it is necessary to circumvent the fact that the sample weights of Y and X ( R 1 and R 2 , which are the row sums of Y and X, respectively) are imposed and are not similar (Supporting Information, Eq. 6 and Eq. 7). This does not meet a crucial constraint in estimating the co-structure between two tables, namely that the samples must be weighted in the same way for the two tables [START_REF] Dray | Co-inertia analysis and the linking of ecological data tables[END_REF]. Hence, an additional common sample weights matrix ( R 0 ) is defined to replace R 1 and R 2 so that the weighted averaging properties of CA can be retained in CoCA [START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF]. In pCoCA, R 0 is equal to the sample weights of the community response table (i.e., R 0 =R 1 ), whereas in sCoCA, Supporting Information,Eq. 10 and Eq. 11). Note that in the specific case of compositional data, the row sums are equal to 1 and R 0 =R 1 =R 2 .

R 0 =( R 1 + R 2 )/ 2 (

Ordination of the structure and assemblage of interacting communities

From sCoCA, ordination diagrams can be made in the usual way by jointly plotting the species scores and sample scores (of each community) for the first axes of the analysis [START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF]. For an optimal representation of this association in a biplot, the species scores of each axis must be multiplied by the quarter root of the eigenvalue of the axis [START_REF] Ter Braak | Interpreting canonical correlation analysis through biplots of structural correlations and weights[END_REF]. From pCoCA, the fit of the response community (viruses) to the predictive community (microalgal hosts) as well as the variations in the composition of communities can be displayed in ordination diagrams (biplots;[START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF]. For instance, the joint plot of sample scores of hosts (table X) with species scores of viruses (table Y) displays the fit of the virus OTUs from the host communities. The simultaneous plot of sample scores of microalgae with the loadings (i.e., the coefficient or importance of variables on the first components of the PLS) of predictor species (i.e., hosts) allows representing microalgae communities and their OTUs. Both types of OTUs can also be jointly displayed using scores of the response species (i.e., viruses) with the loadings of predictor species (i.e., microalgae).

For all these diagrams, axes are optimized to maximize covariance between assemblages rather than to depict associations within individual species matrices, but interpretation can also be carried out as in correspondence analysis according to the barycentre principle, where OTUs are placed at the barycentre (weighted average) of the sample points and, by symmetry, samples at the barycentre (weighted average) of the OTUs (Supporting Information, Eq. 13-16).

Real data application

We illustrate the use of CoCA to study the congruence between microbial communities from NGS data, employing data of [START_REF] Simon | Marked seasonality and high spatial variability of protist communities in shallow freshwater systems[END_REF] on the study of autecology of microbial eukaryotes in shallow freshwater systems, and a data set acquired during our research program (ANR program DECOVIR-12-BSV7-0009) on the monitoring of microalgae and viruses in marine systems.

Case study 1: Microbial eukaryotes

Surface water was sampled monthly from April 2011 to April 2013 in five small shallow freshwater systems [four ponds: Etang des Vallées (EV), La Claye (LC), Mare Gabard (GB), Saint Robert (SR), and one brook: Ru Sainte Anne (RSA)], located at the Natural Regional Park of the Chevreuse Valley (South of Paris, France). These systems were characterized by different local environmental conditions. Briefly, raw genomic sequences were obtained from 18S rDNA fragments, encompassing the V4 hypervariable region, applying 454 pyrosequencing and filtered to remove potential spurious sequences using a local pipeline [START_REF] Simon | Complex communities of small protists and unexpected occurrence of typical marine lineages in shallow freshwater systems[END_REF]. Sequences from all samples were then processed together and clustered into OTUs at 0.98 similarity cutoff using CD hit [START_REF] Fu | CD-HIT: accelerated for clustering the nextgeneration sequencing data[END_REF], and singletons were eliminated, before to assign OTUs to taxonomic groups based on sequence similarity to the PR2 database (Guillou et al., 2013). From this overall OTU table, we used the method applied by [START_REF] Simon | Marked seasonality and high spatial variability of protist communities in shallow freshwater systems[END_REF] to select the most abundant OTUs. Finally, we subdivided the data set by grouping OTUs in two functional groups (autotrophic and heterotrophic) trained on literature information [START_REF] Simon | Marked seasonality and high spatial variability of protist communities in shallow freshwater systems[END_REF][START_REF] Genitsaris | Small-scale variability of protistan planktonic communities relative to environmental pressures and biotic interactions at two adjacent coastal stations[END_REF] Christaki, 2016 and references within), and we obtained one table for autotrophic microbial eukaryotes (n=108, p=122) and another table for microbial heterotrophic eukaryotes (n=108, q=104).

Case study 2: Microalgae-virus system

The data set coming from our research program contains a NGS-based eukaryotic microalgae community table (photosynthetic picoeukaryotes in the class [START_REF] Chen | Amplification of DNA polymerase gene fragments from viruses infecting microalgae[END_REF][START_REF] Clerissi | Unveiling of the diversity of Prasinoviruses (Phycodnaviridae) in marine samples by using high-throughput sequencing analyses of PCR-amplified DNA polymerase and major capsid protein genes[END_REF]. For Mamiellophyceae, the sequence of the V9 region of the 18S rDNA was amplified using primers defined by Amaral-Zettler, McCliment, [START_REF] Amaral-Zettler | A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA genes[END_REF]. The genomic sequences of PolB and V9 region were amplified and sequenced using an Illumina MiSeq platform (GeT-PlaGe, INRA, Castanet-Tolosan, France). Sequences were processed and clustered into OTUs at 0.99 similarity cutoff and singletons were removed using MOTHUR v 1.35.1 [START_REF] Schloss | Introducing mothur: Open-source, platform-independent, communitysupported software for describing and comparing microbial communities[END_REF] for Mamiellophyceae and USEARCH v7 [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF] combined with MUSCLE software [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF] for Prasinovirus. Sequences were then compared against the PR2 database (Guillou et al., 2013) and the NCBI database for Mamiellophyceae and Prasinovirus, respectively, in order to assign OTUs to taxonomic groups based on similarity. Finally, we focused specifically on OTUs assigned to the family Mamiellales of the class Mamiellophyceae, and notably the genus Bathycoccus, Micromonas and Ostreococcus which usually dominate this class in the Gulf of Lion and more generally the picoeukaryotic fraction in other ecosystems [START_REF] Wu | Photosynthetic picoeukaryote assemblages in the South China Sea from the Pearl River estuary to the SEATS station[END_REF][START_REF] Zhu | Mapping of picoeukaryotes in marine ecosystems with quantitative PCR of the 18S rRNA gene[END_REF]. Subsequently, the Prasinovirus data set was limited to OTUs assigned as Bathycoccus viruses (BpVs), Micromonas viruses (MpVs) and Ostreococcus viruses (OtVs). The dominant microalgae and virus OTUs (i.e., ≥ 0.1% of mean relative abundance for at least two samples) were selected to obtained one microalgae table ( n=31, p=67) and one virus table (n=31, q=98).

Statistical analysis

Tables X and Y of autotrophic and heterotrophic microbial eukaryotes, respectively, were subjected to sCoCA. To test the significance of the global co-variation between the two tables, a Monte-Carlo permutation procedure with 9999 permutations was used. In each permutation, sCoCA (by considering all axes) is reapplied to obtain a value of the co-variance between table Y and row-permuted X (so that samples are randomized while preserving the relative abundance of individuals). Note that the choice of the table to be reordered is not important here since we used the symmetric form of the CoCA. A null distribution was estimated from co-variance calculated for the permuted data. The observed co-variance is then compared to the distribution obtained under the null hypothesis. The positions of the samples on ordination axis of each table are then correlated to show the overall level of co-variation between them. For the microalgae-virus data set, both tables were subjected to pCoCA with the SIMPLS algorithm. ter [START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF] suggest that the number of axes used to summarize the data can be selected by "leave-one-out" cross-validation procedure to maximize the cross-validatory fit (%) that measures how well the table X (microalgae in our case) predicts the response table Y (viruses in our case). Working just with these significant axes provides a measure of association between the tables by removing random noise and keeping only the major dimensions of ecological variability. Finally, we combine CoCA and network analysis so that nodes in co-occurrence networks are reordered according to the species scores from the CoCA, and thus from the co-structure between communities. A novel extension of SPIEC-EASI [START_REF] Tipton | Fungi stabilize connectivity in the lung and skin microbial ecosystems[END_REF]) was used to infer the cross-group co-occurrence networks between two data sets. We used the neighborhood (MB) setting and selected the optimal sparsity parameter based on the Stability Approach to Regularization Selection (StARS) [START_REF] Liu | Stability approach to regularization selection (StARS) for high dimensional graphical models[END_REF]. The StARS variability threshold was set to 0.05 for networks built from the two data sets. All statistical analyses were performed with the R software (R Core Team, 2019) and using the cocorresp package [START_REF] Simpson | Cocorresp: Co-correspondence analysis ordination methods[END_REF] for CoCAs and the SpiecEasi package [START_REF] Kurtz | Sparse and compositionally robust inference of microbial ecological networks[END_REF] for co-occurrence networks. Appendix 1 contains a R script and example data (from the case study 1 and the case study 2) allowing users to reproduce the analysis and apply them on their own data sets.

Results

Case study 1: Microbial eukaryotes

The common variance between the two microbial groups computed from the sCoCA explained, significantly ( ), 13.21% of the total variation of autotrophic microbial eukaryotes and 15.22% in heterotrophic microbial eukaryotes. Of the common variance, 39.51% was accounted by the first three axes of the sCoCA (sCoCA axis 1: 18.95%, sCoCA axis 2: 10.55%, sCoCA axis 3: 10.01%). The first three ordination axes of the autotrophic eukaryotes were highly correlated with the first three ordination axes of the heterotrophic eukaryotes (correlations being 0.95, 0.92, and 0.92), demonstrating a high degree of similarity in change between autotrophic and heterotrophic microbial eukaryote assemblages. Communities of autotrophic and heterotrophic microbial eukaryotes covaried along a brook/pond gradient on the first axis (from left to right), and an inter-pond variability on the second axis (Figure 1). Marked differences in the composition of the two communities are visible in joint plots. In the brook system (i.e., RSA), the heterotrophic microbial eukaryote community is mainly composed of fungi, MAST, Labyrinthulida, and Telonema, whereas in pond systems (i.e., EV, LC, MG, SR) Ciliophora, Biocosoecida, Katablepharida, and Choanoflagellida dominated the communities (Figure 1a). Differences between pond systems are explained with higher relative abundance of Biocosoecida, and Katablepharida in EVs and LCs and higher relative abundance of Ciliophera in SRs. Species scores of autotrophic microbial eukaryotes indicated that the patterns in heterotrophic communities are associated to a structure of the autotrophic community (Figure 1b). In the brook system, the autotrophic microbial eukaryote community exhibits high relative abundances of specific OTUs of Bacillariophyceae, Chrysophyte and Cryptophyta. In pond systems, autotrophic communities made up mainly of other specific OTUs of Chlorophyta, Chrysophyte, and Cryptophyta. It is also worth noting that Dinophyta and Streptophyta were found exclusively in pond systems (in particular in MG and SR, respectively).

The community organization of microbial eukaryotes was highlighted from a cross-group co-occurrence network between autotrophic and heterotrophic individuals.

Among 226 dominant autotrophic and heterotrophic OTUs, 204 displayed 274 associations (Figure 2). From these associations between OTUs in the network, more positive (98.5%) than negative associations were inferred. All negative associations occurred between OTUs assigned to Ciliphora for heterotrophic microorganisms and Cryptophyta for autotrophic microorganisms. No clear association patterns can be identified in network from raw tables of autotrophic and heterotrophic microbial eukaryotes (Figure 2a). When the co-occurrence network is reordered according to species scores on the first axis of sCoCA, two modules can be distinguished (Figure 2b).

The first module (i.e., top right corner) is constituted by OTUs exhibiting the higher relative abundances in brook system (i.e., RSA), while OTUs that compose the second module (i.e., bottom left corner) dominate pond systems (i.e., EV, LC, MG, and SR).

Heterotrophic OTUs exhibited major associations with autotrophic OTUs belonging to the same module, with only 1.5% of associations between OTUs from distinct module. A striking pattern is that Chrysophyte is the autotrophic group that contributes most to associations in the two modules (module 1: 71%, module 2: 41%), whereas for the heterotrophic group it is fungi in the module 1 (47%) and Ciliophora in the module 2 (59%). In pond systems (i.e., module 2), surprisingly, fungi are involved in very few associations (7%).

Case study 2: Microalgae-virus system

The cross-validation procedure identified the best pCoCA model based on the first two significant axes (pCoCA axis 1: p=0.001, pCoCA axis 2: p=0.001), in which Mamiellophyceae community predicted 32.02% of the variation in Prasinovirus community. The first two axes accounted for 37.47% (24.26% and 13.21% for axis 1 and 2 respectively) and 44.82% (28.01% and 16.81%) of the variation in the structure of Mamiellophyceae and Prasinovirus community respectively. The biplots indicated that the two communities covaried along a lagoon (samples LAs)/open-sea gradient (samples MAs) on the first axis (from left to right), while a temporal gradient for site LA (intra-site variability) could be identified along the second axis (Figure 3). OtVs have a higher prevalence in the lagoon samples (especially LAs) and coastal samples (SAs) in which Ostreococcus exhibited a high density (Figure 3a Based on the cross-group co-occurrence networks analysis, 67 associations were identified between the major 67 OTUs assigned to one of the three groups of Mamiellophyceae (i.e. Bathycoccus, Micromonas, and Ostreococcus) and the major 98 OTUs assigned to Prasinovirus (i.e. BpVs, MpVs, and OtVs) (Figure 4a). Reordering the co-occurrence network according to the species scores on the first axis of pCoCA highlighted a structure in the network (Figure 4b). The network topology suggests that the identity of OTUs contained in co-occurring groups of viruses and microalgae are related to their respective prevalence along the lagoon/open-sea gradient. Virus OTUs mostly present in lagoon samples have significant associations primarily with microalgae OTUs displaying the higher prevalence in lagoon samples (top right corner, Figure 4b). Similarly, virus OTUs dominating the open-sea samples were mainly associated with microalgae OTUs from open-sea samples (bottom left corner, Figure 4b).

Among the associations contained in the network, 50 were identified between OTUs belonging to an expected host-virus system (i.e. associations Bathycoccus/BpV, Micromonas/MpV, and Ostreococcus/OtV) while the other 17 significant associations were found between OTUs belonging to different host-virus systems. In average 50%, 70.6% and 91.3% of associations found for OTUs of BpVs, MpVs, and OtVs respectively were with OTUs assigned to their respective host group. Within the associations, some single Mamiellophyceae OTUs were associated with many Prasinovirus OTUs and reciprocally. On the other hand, dyads were identified, that is specific associations, in Bathycoccus/BpV, Micromonas/MpV and Ostreococcus/OtV systems. Few negative associations inferred from the observation that those OTUs do not co-occur were found in network (Figure 4b). Interestingly, eight negative associations from a total of nine involved virus OTUs and microalgae OTUs belonging to different host-virus systems.

Discussion

Critical review and guidance papers on the analysis of NGS-based community data [START_REF] Ramette | Multivariate analyses in microbial ecology[END_REF][START_REF] Buttigieg | A guide to statistical analysis in microbial ecology: a community-focused, living review of multivariate data analyses[END_REF][START_REF] Paliy | Application of multivariate statistical techniques in microbial ecology[END_REF] do not mention any direct quantitative method for predicting the composition of one community from another. Co-correspondence analysis [START_REF] Ter Braak | Co-correspondence analysis: a new ordination method to relate two community compositions[END_REF]) fills this gap. At the level of the case study 1, symmetric form of CoCA indicated that heterotrophic microbial eukaryote assemblages in shallow freshwater ecosystems were strongly associated with patterns of autotrophic microbial eukaryotes presence and abundance with links that can be taxon-specific (Figure 1). Our results shown also that the composition of heterotrophic microbial eukaryote community was dominated by fungi in brook system compared to ponds. This is in line with recent observations, based on the estimation of ergosterol level, of a generally higher fungal biomass in river than in ponds [START_REF] Baldy | Microbial dynamics associated with leaves decomposing in the mainstem and floodplain pond of a large river[END_REF]. Higher fungal abundance might potentially be linked to incoming resources from runoff, since in brooks the most important source of imported material is usually deciduous leaves, whose the decomposition processing involved to a large extent fungi [START_REF] Bärlocher | The role of fungi in the nutrition of stream invertebrates[END_REF][START_REF] Webster | Vascular plant breakdown in freshwater ecosystems[END_REF]. To this, the composition of heterotrophic microbial eukaryote community characterizing brook system is associated a specific composition of microbial autotrophs. Such a result suggests that heterotroph community composition exert a control on the composition of autotroph community, and that microbial autotrophs can be driver of microbial heterotrophs.

Regarding the case study 2, the predictive form of CoCA points out different distribution patterns among the three groups of Prasinovirus along the lagoon/open-sea gradient (Figure 3). Note the importance of the dimension reduction step in pCoCA that allows focusing on ecological structures depicted on a limited number of axes and removes random variation from the data. The patterns in Prasinovirus assemblages, with a dominance of OtV in lagoon and coastal samples compared to offshore locations, the inverse distribution for BpV, and MpV exhibiting a wider spatial distribution, are in part a consequence of the presence of their respective hosts in lagoon, coastal and open-sea samples. Indeed, Ostreococcus is known to be abundant in lagoons [START_REF] Subirana | Morphology, genome plasticity, and phylogeny in the genus Ostreococcus reveal a cryptic species, O. mediterraneus sp. nov. (Mamiellales, Mamiellophyceae)[END_REF]), more eutrophic system, compared to Bathycoccus, which is found mainly in oligotrophic areas [START_REF] Vaulot | Metagenomes of the Picolaga Bathycoccus from the Chile coastal upwelling[END_REF][START_REF] Wu | Photosynthetic picoeukaryote assemblages in the South China Sea from the Pearl River estuary to the SEATS station[END_REF]) such as offshore sites (i.e., MA). Micromonas is ubiquitous and particularly present in nutrient-rich environments [START_REF] Not | A single species, Micromonas pusilla (Prasinophyceae), dominates the eukaryotic picoplankton in the western English Channel[END_REF][START_REF] Viprey | Wide genetic diversity of picoplanktonic green algae (Chloroplastida) in the Mediterranean Sea uncovered by a phylum-biased PCR approach[END_REF]. Our findings confirm also the data of [START_REF] Bellec | Abundance, spatial distribution and genetic diversity of Ostreococcus tauri viruses in two different environments[END_REF] showing that OtV are more abundant in lagoon than in the open sea.

Cross-taxon congruence description and evaluation [START_REF] Virtanen | Community concordance between bryophyte and insect assemblages in boreal springs: A broad-scale study in isolated habitats[END_REF][START_REF] Gioria | Evaluation an interpretating vorss-taxon congruence: Potential pitfalls and solutions[END_REF] provides a more comprehensive picture of the community similarity than the richness metrics conventionally used [START_REF] Wolters | Relationship among the species richness of different taxa[END_REF][START_REF] Westgate | Global meta-analysis reveals low consistency of biodiversity congruence relationships[END_REF]. Our results reinforce the need to use CoCA to study the cross-taxon congruence in microbial communities from NGS data. This is all the more important because the high co-correspondence between the two functional groups in microbial eukaryote community may be especially informative given the key ecological role of microbial eukaryotes [START_REF] Caron | Marine protistan diversity[END_REF]. In addition, the study of cross-taxon congruence between Mamiellophyceae and Prasinovirus is of definite interest in marine ecosystems, warmed by climate change, where the expected gradual shift towards small primary producers could render the role of small eukaryotes more important than they are today [START_REF] Morán | Increasing importance of small phytoplankton in a warmer ocean[END_REF]. Microbial eukaryotes are recognized as a significant contributor across various geographical locations of picophytoplankton [START_REF] Worden | Assessing the dynamics and ecology of marine picophytoplankton: the importance of the eukaryotic component[END_REF][START_REF] Jardillier | Significant CO 2 fixation by small prymnesiophytes in the subtropical and tropical northeast Atlantic Ocean[END_REF], which accounts for > 50% of phytoplankton biomass and productivity in marine ecosystems [START_REF] Maranon | Patterns of phytoplankton size structure and productivity in contrasting open-ocean environments[END_REF][START_REF] Teira | Variability of chlorophyll and primary production in the Eastern North Atlantic Subtropical Gyre: potential factors affecting phytoplankton activity[END_REF].

Previous studies have suggested that biotic interactions are the most likely mechanisms underlying cross-taxon congruence at local scales [START_REF] Jackson | Fish and benthic invertebrates: community concordance and community-environment relationships[END_REF][START_REF] Johnson | Spatial congruency of benthic diatom, invertebrate, macrophyte, and fish assemblages in European streams[END_REF], although concordance is also expected from similar responses to environmental gradient [START_REF] Bini | Concordance of species patterns among micro-crustaceans, rotifers and testate amoebae in a shallow pond[END_REF][START_REF] Rooney | Community congruence of plants, invertebrates and birds in natural and constructed shallow open-water wetlands: Do we need to monitor multiple assemblages?[END_REF]. As an important implication, the level of congruence can inform about the structural pattern among interacting groups ( zkan et al., 2014). That being said, Ӧ reordering co-occurrence networks, from the species scores on the first axis of CoCA, allows to identify structural patterns in co-occurrence networks of microbial eukaryote community (Figure 2b) and in the Mamiellophyceae/Prasinovirus system (Figure 4b).

For example, the structure of the microbial eukaryote network is characterized by two modules underlying a brook/pond gradient in the composition of heterotrophic and autotrophic microbial eukaryote assemblages. These differences in OTU composition within the two modules suggest that the food-web structure is different between lotic and lentic ecosystems, and reinforce the differences previously observed at the bacterioplankton level [START_REF] Portillo | Temporal variability in the diversity and composition of stream bacterioplankton communities[END_REF]. A substantial effort has been made in the development of metrics to estimate and test the level of nestedness (e.g., Rodriguez-Guirones & Santamaria, 2006;[START_REF] Ulrich | Null model analysis of species nestedness patterns[END_REF] and modularity (e.g., [START_REF] Banerjee | Netowrk analysis reveals functional redundancy and keystone taxa amongst[END_REF]Dorman & Strauss, 2014) within interaction networks. The order of individuals of two groups in bipartite matrices affects the magnitude of metrics that represent deviations from an idealized state (e.g., perfect nestedness or modularity) [START_REF] Almeida-Neto | A consistent metric for nestedness analysis in ecological systems: reconciling concept and measurement[END_REF]. It has then been advocated that prior to analysis of structure of networks, original bipartite matrices should be reordered to maximize the coherence of individual distributions in rows and columns, so that individuals with most similar links are close together [START_REF] Borgatti | Network analysis of 2-mode data[END_REF][START_REF] Leibold | Coherence, species turnover, and boundary clumping: elements of meta-community structure[END_REF]. These findings taken together with our results validate our proposition to combine CoCA with network analysis to study structural patterns of microbial networks. In addition, in reverse to the expected view of nestedness structure of phage-bacteria network [START_REF] Flores | Statistical structure of host-phage interactions[END_REF], the modular structure of Mamiellophyceae/Prasinovirus network observed in the field underpinned the modularity patterns previously observed in phage-bacteria network from cross-infection experiments [START_REF] Flores | Multi-scale structure and geographic drivers of cross-infection within marine bacteria and phages[END_REF]. Given that the structure of interaction networks is constraint by the coevolutionary processes between species [START_REF] Peralta | Merging evolutionary history into species interaction networks[END_REF], this would lead to account for phylogenetic signals within co-occurrence network [START_REF] Derocles | Biomonitroing for the 21 st century: Integrating next-generation sequencing into ecological network analysis[END_REF]. In this context, it would be possible to disentangle the confounding effect of phylogeny from true biotic interactions by developing a partial analysis (ter Braak, Šmillauer & Dray, 2018) in the context of CoCA to partial out the phylogenetic effect and focus on patterns of co-occurrence that are not related to phylogenetic signal.

All computational methods used to infer networks from NGS data sets produce species co-occurrence networks, where a link between two species represents a significant statistical association (positive or negative) between their abundance. This raises a critical issue about the interpretation of inferred associations [START_REF] Derocles | Biomonitroing for the 21 st century: Integrating next-generation sequencing into ecological network analysis[END_REF], because co-occurrence networks differ from interaction networks constructed on observations of both the species and their interactions [START_REF] Ings | Ecological network -beyond food webs[END_REF]. For instance, all inferred associations between Mamiellophyceae and Prasinovirus belonging to expected host-virus system were positive. These results are consistent with a previous work showing that when parasitism is captured as a significant link in co-occurrence network, it is retrieved as a positive link despite the detrimental effect of parasite on its host [START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF]. This might be explained because the copresence of the host species and the parasite species is necessary for the interaction to occur. Another surprising result is that all negative associations (expected one) were between Mamiellophyceae and Prasinovirus belonging to different host-virus systems. Such negative associations may account for opposite abiotic requirements, since in our case, OTUs of concerned viruses and microalgae had inverse spatio-temporal dynamics.

Positive associations were also found between individuals of different host-virus systems, which could be explained by the increase of a Prasinovirus population that removes, by infection, a major competitor of a co-occurring Mamiellophyceae host of another group. It is important to keep in mind that, although these associations between microalgae and Prasinovirus suggest that they interact, they do not necessarily mean that the co-occurring Mamiellophyceae are the virus hosts, even if they belong to the expected group of hosts (e.g., Bathycoccus, Micromonas or Ostreococcus). Our approach (CoCA combined with network analysis, or other method to infer associations) could be combined with other approaches, e.g. single cell genomics [START_REF] Kalisky | Genomics analysis et the single-cell level[END_REF][START_REF] Martinez-Garcia | Unveiling in situ interactions between marine protists and bacteria through single cell sequencing[END_REF], Epic-PCR [START_REF] Spencer | Massively parallel sequencing of single cells by epicPCR links functional genes with phylogenetic markers[END_REF], to validate the predicted associations in interactions. The justification of the association sign between the two functional groups making up the microbial eukaryote community is also not straightforward, although they can be triggered by ecological interactions or by species abiotic requirements [START_REF] Derocles | Biomonitroing for the 21 st century: Integrating next-generation sequencing into ecological network analysis[END_REF].

In conclusion, the successful application of co-correspondence analysis over two real data sets of microbial communities exhibiting various types of interactions reinforces that resorting to this method for study the distributions, assemblages and interactions between two microbial communities constitutes a highly valuable approach to understand the cross-taxon congruence between microorganisms. A useful consequence of cross-taxon congruence is that the distribution of well-known taxa may provide insight into the processes structuring the distribution of other taxa (e.g., [START_REF] Bilton | How wide to cast the net? Cross-taxon congruence of species richness, community similarity and indicator taxa ponds[END_REF][START_REF] Santi | Simple to sample: vascular plants as surrogate group in a nature reserve[END_REF]. This approach could be used to enhance our understanding of a major problem, the effect of phytoplankton bloom (in particular toxic groups such as cyanobacteria) on the microbial communities and in turn on the ecosystem functions (e.g., [START_REF] Yang | A comprehensive insight into functional profiles of free-living microbial community responses to toxic Akshiwo sanguinea bloom[END_REF][START_REF] Xue | Distinct patterns and processes of abundant and rare eukaryotic plankton communities following a reservoir cyanobacterial bloom[END_REF]. Our findings also demonstrate that the reordering of co-occurrence networks, according to the congruence information extracted from CoCA, allows to investigate the presence of ecological signals in networks. The advantage of this approach is that the complexity of the network is considerably reduced by the non-random placement of nodes in the space in such a way as to improve the aesthetic quality of the representation and consequently its readability, as proposed in good practice of data visualization [START_REF] Spiegelhalter | Visualizing uncertainty about the future[END_REF][START_REF] Kjaergaard | Data visualization: mapping the topical space[END_REF][START_REF] Pocock | The visualization of ecological networks, and their use as a tool for engagement, advocacy and management[END_REF].

Interestingly, the applicability of our approach goes beyond the particular case of data sets with row-sum normalisation (i.e., compositional data). Indeed, CoCA was originally designed to analyze abundance data and is thus able to deal with counts without the need to rarefy data, in accordance with the recent advice against rarefaction [START_REF] Mcmurdie | Waste not, want not: why rarefying microbiome data is inadmissible[END_REF]. It paves the way for further studies to examine the cross-taxon congruence and structural pattern of co-occurrence networks in microbial communities and in turn their effects on the ecosystem functioning.

  Mamiellophyceae) and a NGS-based virus community table (viruses infecting this classof eukaryotic phytoplankton and belonging to the genera Prasinovirus of the family Phycodnaviridae). The data have been acquired across four sites located in northwest Mediterranean Sea (Gulf of Lion) and sampled monthly from March 2013 to April 2014. The Gulf of Lion is characterized by contrasted environments, including eutrophic lagoons connected to the sea, nutrient-rich coastal sites, and oligotrophic open-sea locations. The sample locations included two sites in Leucate lagoon (one coastal site (LB) and another site (LA) at the level of the Grau, i.e. the connection with the sea), a coastal site (SA, marine station included in the French marine monitoring network SOMLIT) and an open-sea site (MA, marine station included in the monitoring network MOOSE). The characterization of Prasinovirus was based on analyzing the partial sequence of the DNA polymerase gene (PolB) amplified using two primer sets

  , b). Conversely, open-sea samples (MAs) were dominated by Bathycoccus which supported virus assemblages dominated by BpVs. Micromonas showed a wider distribution, with a relative contribution of its OTUs both in lagoon samples, coastal samples and open-sea samples, associated with a similar repartition of MpVs (Figure 3a, b).

FiguresFigure 1 .

 1 Figures Figure 1. Ordination biplots of the case study 1, representingof the communities of the positions of sites (open diamond) and species (solid triangle) on the axis 1 × axis 2 factorial plan of the symmetric co-correspondence analysis. (a) Biplot for the heterotrophic microbial eukaryotes and (b) biplot for the autotrophic microbial eukaryotes OTUs obtained from symmetric co-correspondence analysis. EV: Etang des Vallées, LC: La Claye, GB: Mare Gabard, SR: Saint Robert, RSA: Ru Sainte Anne. Heterotrophic and autotrophic microbial eukaryotes OTUs are colored according to the phylogenetic group they belong to. "d" indicates the mesh of the grid.

Figure 2 .

 2 Figure 2. Heatmap of the case study 1, representing the association network between theof heterotrophic and autotrophic microbial eukaryotes (a) before (a) and (b) after (b) the reordering the position of each species in the network according to their network from the species scores on the first axis of symmetric co-correspondence analysis. Bar plots of the relative abundance of (c) autotrophic microbial eukaryotes and (d) heteretrophic microbial eukaryotes. Each OTU is represented by a vertical line partitioned into segments corresponding to its relative abundance in one of five sites. EV: Etang des Vallées, LC: La Claye, GB: Mare Gabard, SR: Saint Robert, RSA: Ru Sainte Anne.

Figure 3 .

 3 Figure 3. Ordination biplots of, concerning the case study 2, ofrepresenting the positions of sites (open diamond) and species (solid triangle) on the axis 1 × axis 2 factorial plan of the predictive co-correspondence analysis. (a) Biplot for the Mamiellophyceae and (b) biplot for the Prasinovirus OTUs obtained from predictive co-correspondence analysis. LB: coastal site in Leucate lagoon, LA: site at the level of the Grau in Leucate lagoon, SA: coastal site, MA: open-sea site. Microalgae and virus OTUs are colored according to the phylogenetic group they belong to. BpV: Bathycoccus viruses, MpV: Micromonas viruses, OtV: Ostreococcus viruses. "d" indicates the mesh of the grid.

Figure 3 .

 3 Figure 3. Heatmap of heterotrophic and autotrophic microbial eukaryotes before (a) and after (b) the reordering the network from the species scores on the first axis of symmetric co-correspondence analysis. Bar plots of the relative abundance of (c) autotrophic microbial eukaryotes and (d) heteretrophic microbial eukaryotes. Each OTU is represented by a vertical line partitioned into segments corresponding to its relative abundance in one of five sites. EV: Etang des Vallées, LC: La Claye, GB: Mare Gabard, SR: Saint Robert, RSA: Ru Sainte Anne.

Figure 4 .

 4 Figure 4. Heatmap of the case study 2, representing theof microalgae-viruses (Mamiellophycaea/Prasinovirus) association network (a) before (a) and (b) after (b) the reordering the position of each species in the network according to their network from the species scores on the first axis of predictive co-correspondence analysis. Bar plots of the relative abundance of (c) microalgae and (d) viruses. Each OTU is represented by a vertical line partitioned into segments corresponding to its relative abundance in one of four sites. LA: site at the level of the Grau in Leucate lagoon, LB: coastal site in Leucate lagoon, SA: coastal site, MA: open-sea site.

Acknowledgments

We are grateful for financial support from the French National Research Agency (DECOVIR ANR-12-BSV7-0009, coordinator Y.D.). We thank Ludwig Jardillier for give us access to genomic data on autotrophic and heterotrophic microbial eukaryotes.

Data and R scripts to reproduce the different analyses of the case study 1 and the case study 2 Mamiellophyceae/Prasinovirus system are available online in the Appendix 1as additional supporting information.

Supporting information

Additional supporting information may be found online in the Supporting Information section at the end of the article.


Susan J Conway
email: susan.conway@univ-nantes.fr

Jack Wright

Costanza Morino

Left frozen block and right, cone shaped molard

Molards are cones of debris that result from the disaggregation of ice-cemented blocks transported by mass movements. Recently, the origin of molards has been directly linked to permafrost degradation on Earth. We have found candidate molards in crater ejecta on Mercury and Mars, which by analogy we link to volatile loss post-emplacement.

Introduction

Distinctive landforms resulting from the degradation of discontinuous mountain permafrost are rare. In our recent research in Iceland, we have uncovered a landform that fills this gap -so called "molards" (Figure 1). These landforms are cones of debris which result from the disaggregation of ice-cemented blocks transported by mass movements [1,2]. We reported for the first time the direct association between the ice cemented block and the resulting molard landform [3]. Permafrost conditions are a prerequisite for forming the ice cement within the blocks, and its subsequent degradation triggers the mass movement leading to the formation of the final molards. In our paper, we show that molards are recognisable in remote sensing data. In addition, we report that the distribution and morphology of molards give insights into landslide dynamics. This new terrestrial work led us to search for such landforms on other planets where the role of volatiles in landscape evolution is under debate. Volatiles are substances that easily change phase, for this work we refer to ices (H 2 O, CO 2 , N 2 , S, Cl) becoming gas or liquid. Our preliminary work has revealed molardlike landforms on Mercury and on Mars.

Mercury: Caloris Ejecta

Conical mounds are found around the Caloris Basin on Mercury. We mapped their locations and found that their distribution is consistent with them being part of the Caloris ejecta blanket. Their shape and the fact that debris from their flanks buries surrounding materials (Figure 2) suggests they continued to evolve long after the basin was formed. Because their flank-slopes and shapes are similar to molards on Earth we suggest a similar formation via volatile loss. In the last few years Mercury has been revealed to be much more volatile-rich than previously thought [e.g. 4], and our discovery of candidate molards in the ejecta of the Caloris Basin suggests that the crust and mantle has been enriched in volatiles since the basin formed ~3.8 billion years ago.

Mars: Hale Crater

We have found candidate molards in the ejecta flows of the one billion year old Hale Crater, which are similar in morphology and spatial distribution to molards found in the Mount Meager debrisavalanche deposits on Earth [5] (Figures 3 &4). Hale impact has already been proposed to have occurred into ice-rich materials [6] and hence it is not unlikely that the ejecta may also have contained ice. The presence of molards suggests that not all the ice was melted by the impact and that the ice itself was contained within the pores of the soil, rather than massive glacial ice. In future work we aim to constrain both the initial ground ice content at the Hale impact, but also the conditions during ejecta emplacement.

Figure 2 :

 2 Figure 2: Conical knob on Mercury. Yellow arrow indicates crater on Caloris Planitia partially obscured by knob material. A subtle notch in the knob is visible above this crater.

Figure 3 :

 3 Figure 3: Oblique view of candidate molards in Hale crater ejecta. Cone basal diameter is about 150 m.

Figure 4 :

 4 Figure 4: Comparison between candidate molards in Hale crater ejecta, left and molards in Mount Meager, right.

Acknowledgements

SJC is grateful for the financial support of CNES in support of her HiRISE and BepiColombo work. This work acknowledges the financial support from Région Pays de la Loire, project GeoPlaNet (convention N° 2016-10982). JW is funded by Planmap (EU Horizon 2020 grant no. 776276).