Storme Zaviar De Scally

Samuel Chaffron

Thulani Peter Makhalanyane
email: thulani.makhalanyane@up.ac.za

Polar opposites; bacterioplankton susceptibility and mycoplankton resistance to ocean acidification

Keywords: bacterioplankton, mycoplankton, acidification, Southern Ocean, diversity, enzymatic activity, carbon, nitrogen, phosphorus

associated with nitrogen acquisition. This study adds to our understanding of the effects of ocean acidification on microbial communities, particularly within an environment expected to be largely effected by future anthropogenically driven climate change.

ABSTRACT

Microorganisms form the basis of ocean ecosystems yet the effects of perturbations such as decreasing pH on microbial community structure, interactions and functionality remain compared to multicellular organisms. Using an experimental manipulation of Southern Ocean seawater, we subjected bacterioplankton and mycoplankton to artificial pH decreases, which are predicted to occur in the future. We show that acidification led to substantial increases of bacterioplankton diversity, while in contrast it had no effect on mycoplankton diversity. Our analyses revealed a loss of putative keystone taxa and a decrease in predicted community interactions as a response to lower pH levels.

Bacterioplankton shifted from generalist to specialist community members, suggesting a specific stress response to unfavourable conditions. In addition, enzyme activities involved in nitrogen acquisition were lower at reduced pH levels, suggesting altered organic matter cycling in a more acidic ocean. Our findings suggest that bacterioplankton and mycoplankton may respond differentially to future ocean acidification, with potentially negative impacts on community structure and biogeochemical cycling in the Southern Ocean.

IMPORTANCE

Oceans absorb the majority of anthropogenically produced CO 2 , the consequence of which is ocean acidification, a phenomenon already negatively impacting key marine organisms. Marine microbial communities form the basis of ocean food webs by generating nutrients for higher trophic levels, yet the response of these key microbial drivers to acidification remains unclear. This knowledge deficit is particularly true for understudied marine ecosystems such as the Southern Ocean. Using a mesocosm approach, we found that acidification severely impacts microbial community stability, by altering bacterioplankton community structure, reducing network complexity, and augmenting enzyme activities

INTRODUCTION

The ocean absorbs roughly 30% of anthropogenic carbon dioxide emissions, which leads to changes in oceanic carbonate chemistry, resulting in both carbonation and acidification [START_REF] Orr | Anthropogenic ocean acidification over the twenty-first century and its impact on calcifying organisms[END_REF][START_REF] Zeebe | Carbon emissions and acidification[END_REF][START_REF] Doney | Ocean acidification: the other CO2 problem[END_REF][START_REF] Adler | Productivity Is a Poor Predictor of Plant Species Richness[END_REF][START_REF] Kottmeier | Acidification, not carbonation, is the major regulator of carbon fluxes in the coccolithophore Emiliania huxleyi[END_REF]. Globally, ocean pH is predicted to decrease by as much as 0.4 and 0.8 units within the next 100 and 300 years, respectively [START_REF] Caldeira | Oceanography: anthropogenic carbon and ocean pH[END_REF][START_REF] Orr | Anthropogenic ocean acidification over the twenty-first century and its impact on calcifying organisms[END_REF], a change which is nearly eight times higher than what has been observed within the last 25 million years [START_REF] Pearson | Atmospheric carbon dioxide concentrations over the past 60 million years[END_REF]. Ocean acidification (OA) can negatively impact marine macro-organisms, such as calcifiers, adversely affecting their growth and physiology [START_REF] Doney | Ocean acidification: the other CO2 problem[END_REF][START_REF] Gazeau | Impacts of ocean acidification on marine shelled molluscs[END_REF][START_REF] Brandenburg | Meta-analysis reveals enhanced growth of marine harmful algae from temperate regions with warming and elevated CO2 levels[END_REF]. However, the extent to which decreasing pH affects diverse bacterioplankton and mycoplankton communities remains largely unclear.

A null hypothesis proposed by Joint and colleagues (2011b) states that due to the ubiquitous abundance, fast generation time and genetic plasticity of microorganisms, microorganisms will rapidly adapt to future pH perturbations. However, despite the large genetic flexibility microorganisms possess, a physiological tipping point, where individual organisms are no longer able to adapt to perturbations, may be reached (Allison and Martiny, 2008b;[START_REF] Hutchins | Microorganisms and ocean global change[END_REF]. There is considerable uncertainty and disagreement regarding the impact of acidification on bacterioplankton and mycoplankton communities, with reports demonstrating both resistance and susceptibility. Evidence of microbial community resistance to acidification include: unchanged bacterial cell abundance [START_REF] Grossart | Testing the effect of CO 2 concentration on dynamics of marine heterotrophic bacterioplankton[END_REF][START_REF] Allgaier | Coupling of heterotrophic bacteria to phytoplankton bloom development at different pCO 2 levels: a mesocosm study[END_REF][START_REF] Krause | Small changes in pH have direct effects on marine bacterial community composition: a microcosm approach[END_REF][START_REF] Newbold | The response of marine picoplankton to ocean acidification[END_REF], consistent composition [START_REF] Lindh | Consequences of increased temperature and acidification on bacterioplankton community composition during a mesocosm spring bloom in the Baltic Sea[END_REF][START_REF] Roy | Ocean acidification shows negligible impacts on high-latitude bacterial community structure in coastal pelagic mesocosms[END_REF][START_REF] Oliver | Marine bacterial communities are resistant to elevated carbon dioxide levels[END_REF], unaffected or increased heterotrophic production [START_REF] Grossart | Testing the effect of CO 2 concentration on dynamics of marine heterotrophic bacterioplankton[END_REF][START_REF] Arnosti | Dynamics of extracellular enzyme activities in seawater under changed atmospheric pCO2: a mesocosm investigation[END_REF] and increased nitrogen fixation (Joint et al., 2011a;[START_REF] Piontek | Ocean acidification impacts on nitrogen fixation in the coastal western Mediterranean Sea[END_REF]. Susceptibility has been seen through, alterations in bacterial and fungal cell numbers [START_REF] Krause | Marine fungi may benefit from ocean acidification[END_REF][START_REF] Endres | Stimulated bacterial growth under elevated pCO2: results from an off-shore mesocosm study[END_REF][START_REF] Engel | Impact of CO 2 enrichment on organic matter dynamics during nutrient induced coastal phytoplankton blooms[END_REF], compositional shifts [START_REF] Allgaier | Coupling of heterotrophic bacteria to phytoplankton bloom development at different pCO 2 levels: a mesocosm study[END_REF][START_REF] Krause | Small changes in pH have direct effects on marine bacterial community composition: a microcosm approach[END_REF][START_REF] Zhang | Response of bacterioplankton community structure to an artificial gradient of pCO2 in the Arctic Ocean[END_REF] with changes towards pathogenic members [START_REF] Vega-Thurber | Metagenomic analysis of stressed coral holobionts[END_REF][START_REF] Reich | Impacts of a reduction in seawater pH mimicking ocean acidification on the structure and diversity of mycoplankton communities[END_REF]) and decreased nitrification rates [START_REF] Beman | Global declines in oceanic nitrification rates as a consequence of ocean acidification[END_REF]. Moreover, the increased expression of proton pump genes identified in a metatranscriptomic study [START_REF] Bunse | Response of marine bacterioplankton pH homeostasis gene expression to elevated CO2[END_REF] suggests that acidification may be energetically costly for microbial cells, with possible negative implications for microbial functionality and biogeochemical cycling. These results emphasise the need for further investigations into the effects of acidification on complex and diverse bacterioplankton and mycoplankton assemblages [START_REF] Hutchins | Microorganisms and ocean global change[END_REF].

Although the Southern Ocean (SO) covers less than 20% of the global ocean it is a highly productive ecosystem, facilitating the inter-ocean exchange of nutrients which drives over 70% of ocean productivity [START_REF] Mayewski | State of the Antarctic and Southern Ocean climate system[END_REF]. While the SO has naturally low pH, it remains vulnerable to anthropogenic induced acidification due to its' low temperature, the upwelling of CO 2 enriched deep waters, and the fact that approximately 40% of the global oceanic uptake of CO 2 occurs in the SO [START_REF] Riebesell | Lessons learned from ocean acidification research[END_REF]. Several studies have shown that acidification results in increased resilience in diatoms [START_REF] Valenzuela | Ocean acidification conditions increase resilience of marine diatoms[END_REF] and functional biodiversity loss [START_REF] Teixidó | Functional biodiversity loss along natural CO 2 gradients[END_REF]. In contrast to other environments [START_REF] Newbold | The response of marine picoplankton to ocean acidification[END_REF][START_REF] Chen | Shifts in microbial community composition and function in the acidification of a lead/zinc mine tailings[END_REF][START_REF] Huggett | Future warming and acidification result in multiple ecological impacts to a temperate coralline alga[END_REF][START_REF] James | Elevated pCO2 alters marine heterotrophic bacterial community composition and metabolic potential in response to a pulse of phytoplankton organic matter[END_REF], the precise responses of bacterioplankton and mycoplankton remain largely speculative and often unexamined in pristine marine environments. This knowledge is crucial for understanding the effects of climate change in regions most sensitive to future environmental change.

To address this deficit, we assessed the effect of future predicted pH levels (based on IPCC projections for the years 2100 and 2300, respectively) on both SO bacterioplankton and mycoplankton diversity, composition, predicted interactions, and functionality using mesocosm experiments. We predicted that; (1) ocean acidification will decrease both bacterioplankton and mycoplankton diversity and will significantly alter community structure; (2) microbial functionality will be both positively and negatively affected by ocean acidification in terms of carbon and nitrogen/phosphorus decomposition, respectively, and (3) the total number of covariations and subsequently network stability will decrease with lower pH levels.

RESULTS

Microbial diversity and community structure in control and acidified mesocosms

Alpha diversity for bacteria and archaea was significantly higher in the low and medium pH groups in comparison to the control (ANOVA, p-value < 0.001), and both Shannon diversity (Fig. 1a) and OTU richness (Fig. 1b) significantly correlated with pH. On the contrary, fungal alpha diversity was not significantly different between the three treatment groups (ANOVA, p-value > 0.05), and did not correlate significantly with pH (Fig. 1ab). Bacterial and archaeal community structure was significantly affected by pH treatment (PERMANOVA, p-value < 0.05, betadisper p-value > 0.05), and clustered according to treatment type (Fig. 2a). In contrast to bacteria and archaea, fungal community structure was not significantly influenced by pH treatment (PERMANOVA, p-value < 0.05, betadisper pvalue < 0.05), as no distinct clusters were seen according to treatment group (Fig. 2b). Thus, bacterioplankton and mycoplankton community composition clearly display differential responses to acidification. Whereas several bacterioplankton genera either increased or decreased in abundance under acidified conditions, mycoplankton composition remained stable, with the exception of Agariomyctes. Specific genera and species differences are further described in Supporting Information (Fig. S3).

Functional community responses to ocean acidification

Extracellular enzymatic activities were detected in all samples, where the activities of carbon and phosphorus acquiring enzymes did not vary significantly with pH treatment (Fig. 3). In contrast, nitrogen acquiring enzymes NAG and LAP significantly decreased at lower pH. Enzymes involved in carbon acquisition, namely BG and BX, showed no significant increase in response to pH treatment (ANOVA, p-value > 0.05). AP, involved in phosphorus acquisition, showed no significant increase or decrease in activity in response to the pH treatment (ANOVA, p-value > 0.05). NAG activity decreased in the low group across sampling days, although this was not significant from the control group (ANOVA, p-value > 0.05). Notably, LAP was shown to be significantly lower in the medium and low pH treatment groups relative to the control at day 21 of the mesocosm experiment (ANOVA, pvalue < 0.05). Activities were correlated to pH levels using both linear regression (Fig. 3) and Spearman's rho. Both results indicated that only NAG was significantly positively associated with pH (r = 0.42, p-value < 0.05; Fig. 3a), whereby NAG activity decreased at lower pH levels. LAP activity, although significantly different between treatment and control groups, could not be predicted by pH, likely due to the lower activity in the medium pH group than the low pH group.

Community structures and system level changes associated to ocean acidification

Table S2 shows topological features of bacterial and archaeal community networks.

All networks exhibited scale-free characteristics, typical of biological association networks, were small world, modular, and significantly different from generated random networks (see Supporting information for details).

Network structure, complexity and community composition differed markedly between the three treatment groups (Table S2, Fig. 4). The medium network was the largest followed by the low and control networks, showing that acidification resulted in larger networks. Interestingly, although the medium network displayed the highest diversity and connectivity, the low network was the least complex (Table S2). The medium network contained the greatest number of edges, the largest average degree (i.e. average number of edges per node in the network) and the shortest harmonic geodesic distance (i.e. distance between all nodes in the network).

The number of modules in the medium and low networks were greater than the control network (Table S2). Modules were largest and most connected in the medium network, followed by the control and the low network, reflecting patterns observed from overall network structure. Intriguingly, the percentage of positive associations increased in response to acidification, where the low network had 77.82% positive associations, the medium network contained 77.25% and the control network only 74.70%, highlighting that reduced pH levels resulted in a slight increase in positive associations as opposed to coexclusions (Fig. 4, Table S2). The overall taxonomic composition of the largest modules at both phylum and class level did not differ markedly between the three treatment groups, although some modules exhibited a dominance of different taxonomic groups (Fig. 4).

Alphaproteobacteria, Bacteroidetes and Gammaproteobacteria were the dominant taxa in almost all modules and were also the dominant taxa in the overall bacterioplankton community. Actinobacteria occurred more frequently in the medium and low networks, and Euryarchaeota was found exclusively in the control network (Fig. S2 and Table S1).

Unassigned sequences constituted a relatively high proportion in several modules, and were highest in the low network, suggesting that cryptic unknown species may increase in abundance in response to acidification in the SO. Furthermore, no modules consisted of a single class or phylum, showing a low degree of phylogenetic assortativity within the association networks, which suggest a potential high capacity for biotic interactions in SO microbial communities.

Network topological metrics can be computed for each node/OTU by assessing the within-module connectivity (Zi) as well as the between-module connectivity (Pi), which can indicate the potential ecological role of each OTU. These parameters can then be used to classify each OTU into four network-based categories; peripherals, connectors, module hubs and network hubs (Fig. 5a). Ecologically, peripherals may be considered as specialists, while connectors and module hubs may represent generalists and network hubs could be classified as super generalists. Additionally, the last three can be regarded as potential keystone taxa, on the sole basis of their high connectivity, since the overall network structure would be significantly altered if they were removed.

The majority of OTUs in all three networks were classified as peripherals, with most edges occurring within their own modules (Fig. 5a). Several peripheral OTUs had a Pi value of 0, indicating that these nodes were connected solely to other OTUs within the same module (74% of OTUs in the control network, 70% in the medium network and 84% in the low). No network hubs (super generalists) were detected in any network. Interestingly, the relative abundance of several peripheral taxa changed in response to acidification (Fig. 5b, see Supplemental Information for details).

The number and taxonomic affiliation of module hubs and connectors detected (Table S3) differed noticeably between the three networks. Five module hubs were detected in the control network, five in the medium network and only three in the low network (Fig. 4,5a).

All control hubs consisted of phylum Bacteroidetes OTUs and unassigned sequences, with the orders Flavobacteriales and Saprospirales. Medium hubs belonged to classes Flavobacteriia, Gammaproteobacteria and Alphaproteobacteria. Low hubs consisted of Alphaproteobacteria, Flavobacteriia and unassigned bacteria. The relative abundance of module hubs ranged from 0.01% to 0.74%. The control network contained five connectors, the medium network contained three and the low network had only one. Control and medium connectors consisted of Alphaproteobacteria (order Rhodobacterales and Rickettsiales, respectively) and Flavobacteriia (order Flavobacteriales), whereas the low connector belonged to Bacteroidetes (order Cytophagales). Connector relative abundance ranged from 0.005% to 3.24%, with all but one OTU with an abundance of less than 1%, underlying the disconnection between species centrality and relative abundance. As hubs and connectors can be considered potential keystone taxa, Alphaproteobacteria and Flavobacteriia represented prominent putative keystone taxa within all networks, although differences between networks can be detected at finer taxonomic resolution (Table S3).

DISCUSSION

Diversity does not equal stability

Microbial community stability amidst perturbations is thought to be promoted by increased diversity, due to the link between diversity and ecosystem health and multifunctionality [START_REF] Girvan | Bacterial diversity promotes community stability and functional resilience after perturbation[END_REF][START_REF] Delgado-Baquerizo | Microbial diversity drives multifunctionality in terrestrial ecosystems[END_REF][START_REF] Shade | Diversity is the question, not the answer[END_REF].

Similarly, loss of community diversity following environmental disturbance is often linked to community susceptibility. However, it has recently been suggested that assumptions of increased diversity leading to increased stability are false, at least for some ecosystems, as the role of diversity is dependent on the environment or ecosystem and disturbance itself [START_REF] Shade | Diversity is the question, not the answer[END_REF]. Community members are not necessarily interacting mutualistically in diverse compositions, and taxonomic diversity does not absolutely equate to functional diversity [START_REF] Louca | Decoupling function and taxonomy in the global ocean microbiome[END_REF]. In this study, we found that bacterioplankton alpha diversity is significantly higher in acidified compared to ambient mesocosms, in contrast to what has been hypothesised. Simultaneously, bacterioplankton community structure was significantly altered in response to CO 2 , with the relative abundances of several key community members significantly changing. We suggest that increased bacterioplankton diversity is an indication of community susceptibility, or "scramble" rather than stability in the Southern Ocean. The high diversity identified for bacterioplankton communities highlight that polar marine regions, and the SO specifically, contain high levels of unexplored diversity [START_REF] Dickinson | Extremophiles in an Antarctic marine ecosystem[END_REF].

Mycoplankton resistance to acidification

In contrast to our initial hypothesis, we did not observe a significant effect of ocean acidification on mycoplankton diversity. Instead, fungal community diversity and structure was stable, with the exception of Agaricomycetes, which increased in relative abundances at lower pH levels (Table S3). Agaricomycetes is a class of Basidiomycota that contains several yeast species [START_REF] Jones | Fifty years of marine mycology[END_REF][START_REF] Jones | Classification of marine ascomycota, basidiomycota, blastocladiomycota and chytridiomycota[END_REF], suggesting an increase in abundance of adapted fungal morphologies in response to acidification. However, the increase of Agaricomycetes phylotypes did not significantly affect overall fungal diversity or community structure. We argue that mycoplankton may contain a higher genetic plasticity and degree of physiological adaptation than bacterioplankton and are therefore resistant to future acidification. Culture dependant studies have shown that both terrestrial and marine fungi can rapidly adapt to a variety of environmental conditions, such as hot and cold temperatures and anoxia [START_REF] Richards | Marine fungi: their ecology and molecular diversity[END_REF][START_REF] Tsuji | Cold adaptation of fungi obtained from soil and lake sediment in the Skarvsnes ice-free area, Antarctica[END_REF][START_REF] Rédou | Species richness and adaptation of marine fungi from deep-subseafloor sediments[END_REF]. Additionally, fungi have larger pH tolerance ranges and optima for growth than bacteria [START_REF] Rousk | Soil bacterial and fungal communities across a pH gradient in an arable soil[END_REF].

Moreover, fungi possess chitin-rich cell walls [START_REF] Taylor | Multi-year assessment of coastal planktonic fungi reveals environmental drivers of diversity and abundance[END_REF], which may provide enhanced resistance. Source habitat significantly influences mycoplankton community structure [START_REF] Panzer | Identification of habitatspecific biomes of aquatic fungal communities using a comprehensive nearly full-length 18S rRNA dataset enriched with contextual data[END_REF], further suggesting that SO fungal communities may be adapted to the pelagic environment and therefore less effected by pH variations. The lack of a significant mycoplankton community response in SO waters to pH levels predicted in the next 300 years, suggests the resistance of marine fungi to acidification in these systems.

Altered organic matter cycling in a future acidic ocean

Microbially produced extracellular enzymes break down complex organic substrates, mediating heterotrophic carbon and nutrient recycling in marine environments [START_REF] Arnosti | Patterns of microbially driven carbon cycling in the ocean: links between extracellular enzymes and microbial communities[END_REF]. Enzymes measured for nutrient acquisition showed activity in all samples, indicating active bacterioplankton and/or mycoplankton heterotrophic communities within surface SO waters, since members of both assemblages are capable of producing each enzyme [START_REF] Jones | Fifty years of marine mycology[END_REF][START_REF] Arnosti | Patterns of microbially driven carbon cycling in the ocean: links between extracellular enzymes and microbial communities[END_REF][START_REF] Hoarfrost | Heterotrophic extracellular enzymatic activities in the atlantic ocean follow patterns across spatial and depth regimes[END_REF]. The observed results of exoenzyme activities were congruent with our original hypothesis. Carbon and phosphorus degrading enzymes showed no substantial shifts in activity in response to pH treatment, suggesting constant, and not elevated carbon and phosphorus utilisation. Nitrogen acquiring enzymes, LAP and NAG, however, were negatively associated with pH, implying altered heterotrophic nitrogen recycling in a future more acidic ocean.

LAP and NAG activity responses contrasted to previous studies, which have found LAP activity to increase [START_REF] Piontek | Ocean acidification impacts on nitrogen fixation in the coastal western Mediterranean Sea[END_REF][START_REF] Endres | Stimulated bacterial growth under elevated pCO2: results from an off-shore mesocosm study[END_REF] or remain unchanged [START_REF] Yamada | Impact on bacterial activities of ocean sequestration of carbon dioxide into bathypelagic layers[END_REF][START_REF] Burrell | Assessing approaches to determine the effect of ocean acidification on bacterial processes[END_REF]. The decreased LAP and NAG activities suggest that nitrogen acquisition by bacterioplankton may be negatively impacted by pH, possibly due to increased physiological stress placed on microbial cells in a more acidic ocean. [START_REF] Bunse | Response of marine bacterioplankton pH homeostasis gene expression to elevated CO2[END_REF] found that bacterioplankton respond to a pH decrease of 0.2 units, applied over a period of only nine days, by increasing the expression of pH homeostasis genes. The exposure to lower extracellular pH likely increased physiological stress, resulting in intensified energy expenditure on cell maintenance, at the expense of bacterial growth [START_REF] Bunse | Response of marine bacterioplankton pH homeostasis gene expression to elevated CO2[END_REF]. Therefore, it is likely that the microbial communities producing LAP and NAG were negatively impacted by acidification. Reduced pH may have caused a loss of community members responsible for producing the nitrogen acquiring enzymes (Allison and Martiny, 2008a), ultimately decreasing cellular growth and fitness. Ultimately, microbial communities are essential produces of enzymes responsible for recycling organic matter, especially within nutrient deficient regions such as the SO. Although further research is necessary, alterations in bacterioplankton communities and their functions within the oligotrophic SO suggest negative future effects on ocean nutrient cycling.

Bacterioplankton community network shifts as a response to acidification

Network topology, size and structure shifted markedly in response to reduced pH, with prominent differences between the control, medium and low bacterioplankton networks, congruent with our third hypothesis. Altered network structure and the loss of potential keystone taxa at lower pH levels suggest that marine microbial communities will be impacted by future acidification. Modules within networks may represent specific ecological niches [START_REF] Olesen | The modularity of pollination networks[END_REF][START_REF] Dupont | Stability of modular structure in temporal cumulative plantflower-visitor networks[END_REF], and high interconnectedness between modules may indicate communication or co-operation between niches. The control group, although containing a large number of modules, did not contain a high number of connections between modules, and several between module connections were negative associations (Fig. 4). The large number of modules observed in the control group suggests that SO bacterioplankton communities form distinct niches and that both intermodular and intramodular predicted interactions may represent a balance between competition and mutualism. For example, microorganisms compete for nutrients and exhibit different growth constraints [START_REF] Steele | Marine bacterial, archaeal and protistan association networks reveal ecological linkages[END_REF][START_REF] Eiler | Coherent dynamics and association networks among lake bacterioplankton taxa[END_REF][START_REF] Brum | Seasonal time bombs: dominant temperate viruses affect Southern Ocean microbial dynamics[END_REF] and also communicate via mechanisms such as quorum sensing [START_REF] Hmelo | Quorum sensing in marine microbial environments[END_REF]. This study supports the notion that SO bacterioplankton communities are highly complex with a multitude of putative biotic interactions.

The higher network connectivity and complexity in the medium pH group suggests that bacterioplankton responded to a pH decrease of 0.4 units, associated with an increase in network stability (Fig. 4). Indeed, pH has been shown to alter community network structure [START_REF] Barberan | Using network analysis to explore cooccurrence patterns in soil microbial communities[END_REF], resulting in larger and more complex phylogenetic and functional networks in response to elevated soil CO 2 levels [START_REF] Zhou | Functional molecular ecological networks[END_REF][START_REF] Zhou | Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2[END_REF].

Furthermore, network stability is thought to be promoted by increased modularity, as disturbances may spread more slowly through a modular than a non-modular network [START_REF] Olesen | The modularity of pollination networks[END_REF]. Module 4 contained the highest number of intra-module connections, where the majority of co-occurrences were positive (Fig. 4). The high intra-module connectivity is likely indicative of cooperation or mutualism between phylogenetically distant species within a specific niche. In contrast, several intermodular connections were negative, suggesting that increased physiological stress may result in more competitive interactions between niches at lower pH levels within the SO. This competition may be due to increased physiological stress on weaker community members, leading to the development of specialists with more efficient mechanisms for using the available nutrients. However, more research is needed to determine the cause of competition observed. Ultimately, future pH decreases of 0.4 units may promote increased bacterioplankton community interactions and network stability.

In contrast, the low network was the least complex, with reduced intermodular connectivity. This may suggest that microbial niches or guilds decreased the amount of interactions compared to the control and medium networks. This decreased cooperation may be due to the increased physiological stress placed on bacterioplankton communities at pH levels expected in the next 300 years. Acidification can promote direct physiological changes of microbial communities, whereby the increased energetic cost associated with cellular maintenance may lead to reduced microbial growth and fitness [START_REF] Gilbert | Detection of large numbers of novel sequences in the metatranscriptomes of complex marine microbial communities[END_REF][START_REF] Bunse | Response of marine bacterioplankton pH homeostasis gene expression to elevated CO2[END_REF]. Decreased cellular fitness may result in altered functionality [START_REF] Bunse | Response of marine bacterioplankton pH homeostasis gene expression to elevated CO2[END_REF] eventually perturbing community interactions. The increased stress induced by a pH decrease of 0.8 units likely resulted in the substantial alterations of the network structure, with potentially deleterious effects on bacterioplankton community structure in a future acidified ocean scenario.

Keystone bacterioplankton in the SO and responses to acidification

Several potential keystone taxa were identified in all three networks, defined by their integral topological role in maintaining network structure and stability [START_REF] Deng | Molecular ecological network analyses[END_REF].

The number of connectors decreased in the medium network, suggesting a switch between single microorganisms responsible for connecting niches to several community members linking different modules. The number of module hubs and connectors decreased dramatically in the low network. Loss of keystone members in a network decreases network stability and can lead to a collapse of network structure due to their high level of connectivity [START_REF] Dunne | Food-web structure and network theory: the role of connectance and size[END_REF][START_REF] Olesen | The modularity of pollination networks[END_REF]. The observed decrease in keystone nodes in the low network therefore correlates to the less stable network structure. Moreover, the reduced frequency of module hubs and connectors yet again suggests a lower frequency of community interactions or relationships in response to extreme levels of acidification.

The constructed networks are based on co-abundances and do not necessarily represent true biological interactions [START_REF] Shi | The interconnected rhizosphere: high network complexity dominates rhizosphere assemblages[END_REF][START_REF] Röttjers | Can we predict keystones?[END_REF]. However, co-occurrence networks have been instrumental for investigating the structure of microbial assemblages across biomes [START_REF] De Menezes | Network analysis reveals that bacteria and fungi form modules that correlate independently with soil parameters[END_REF][START_REF] Murdock | Hydrothermal vent protistan distribution along the Mariana arc suggests vent endemics may be rare and novel[END_REF], in the detection of true viral host interactions [START_REF] Chow | Top-down controls on bacterial community structure: microbial network analysis of bacteria, T4-like viruses and protists[END_REF][START_REF] Soffer | Phage-bacteria network analysis and its implication for the understanding of coral disease[END_REF], the association of phylotype shifts to specific environmental parameters and functions [START_REF] Valverde | Cyanobacteria drive community composition and functionality in rock-soil interface communities[END_REF][START_REF] Van Goethem | Cyanobacteria and Alphaproteobacteria May Facilitate Cooperative Interactions in Niche Communities[END_REF], as well as in the detection of metabolic dependencies resulting in the culture of previously uncultivable microorganisms [START_REF] Duran-Pinedo | Correlation network analysis applied to complex biofilm communities[END_REF].

Furthermore, networks provide valuable insight into complex community structures and stability, particularly in the face of perturbations, and give information on community organisation not portrayed by diversity metrics alone [START_REF] Zhou | Functional molecular ecological networks[END_REF]. Thus, our results provide a framework for complex community predicted interactions in response to acidification (Fig. 6) where individual species interactions may be experimentally tested in future studies.

In summary, we have shown that ocean acidification significantly alters bacterioplankton, but not mycoplankton, communities. Exoenzyme activities involved in nitrogen acquisition were substantially reduced at lower pH levels, suggesting a diminished availability of nutrients in surface waters and amended nutrient recycling. Additionally, bacterioplankton network structure and keystone taxa were significantly impacted by increased CO 2 . We propose that future acidification may have adverse consequences on microbial communities and biogeochemical cycling, particularly within regions most sensitive to future climatic changes.

MATERIALS AND METHODS

Sampling and mesocosm experiment

Seawater was collected during the Marion Island Relief Cruise aboard the RV SA Agulhas II from April 2016 to May 2016. Approximately 500 L of seawater was collected from surface waters (5 m) using the underway system at a single site (46°00.065' S, 38°00.37' E). The water was distributed into nine mesocosms (50 L volume) with the sterile drums filled to near capacity to keep the headspace small (volume 5 -10ml) in order to minimize the exchange of CO 2 . The mesocosms consisted of three biological replicates for each pH treatment group, where treatments are defined as follows: a control (i.e. pH 8.0) which is the natural pH of the environment, a medium pH, which mimics predicted pH decreases within the next 100 years (0.4 units, i.e. pH 7.6), as well as a low pH treatment which exceeds the predicted decrease in pH within the next millennium (i.e. pH 7.2). All values were allowed to deviate by ± 0.05 pH units in order to account for natural variability within the SO environment [START_REF] Kapsenberg | Near-shore Antarctic pH variability has implications for the design of ocean acidification experiments[END_REF]. The mesocosms were exposed to light on a 12 hr day/night cycle, where light was allowed to penetrate the top of the mesocosms in order to mimic the natural light conditions at 5 m. Mixing of mesocosms was performed via consistent O 2 bubbling (HAILEA Air Compressor, Guangdong, China), as well as manual stirring throughout the experiment. Each treated mesocosm contained a pH probe (MODEL PH-301C, SAGA, Tainan City, Taiwan, pH accuracy ± 0.1% + 2 digits) which consistently recorded the seawater pH. Physico-chemical variables within each mesocosm such as temperature, salinity, conductivity and pH, were also manually recorded twice a day throughout the duration of the experiment using a HANNA probe (MODEL HI 98195, Padova, Italy, pH accuracy ± 0.02 units, temperature accuracy ± 0.15°C, conductivity and salinity accuracy ± 1% of reading). Changes in pH were induced via CO 2 bubbling and the experiment was run for a total of 21 days, during which pH changes were induced on the first day (day 0) of the experiment. Both O 2 and CO 2 bubbling were introduced into mesocosms via sterile plastic tubing, where the amount of gas was controlled using a regulator. Samples were collected every 7 days at noon, including at the start (day 0) of the experiment. Physicochemical variables such as temperature, salinity, and nutrients were recorded throughout the experiment (see Fig S1).

Enzymatic activities

Seawater collected from each sampling point in the mesocosm experiment was assayed for the activity of 5 enzymes, including β-glucosidase (BG), β-xylosidase (BX), β-Nacetyl-glucosaminidase (NAG), alkaline phosphatase (AP), and leucine aminopeptidase (LAP). The potential activities of these enzymes were tested using fluorescent assays, and enzyme activities were performed and calculated as described previously [START_REF] Sinsabaugh | Enzymic analysis of microbial pattern and process[END_REF][START_REF] Papanikolaou | Nitrogen deposition, vegetation burning and climate warming act independently on microbial community structure and enzyme activity associated with decomposing litter in low-alpine heath[END_REF][START_REF] Hoarfrost | Heterotrophic extracellular enzymatic activities in the atlantic ocean follow patterns across spatial and depth regimes[END_REF]. Briefly, 200 µL of seawater was added to 96-well black flat-bottom microplates (Greiner Bio One, Frickenhausen, Germany), where four replicate wells were used per sample. A concentration of between 1 mM to 1.5 mM of each enzyme substrate was then added to the 96 well plates. Plates were incubated for 2 hours at 4°C in the dark. Fluorescence was measured from the black 96 well plates using a Spectramax® Paradigm Multi-Mode Microplate Reader (Molecular Devices, USA) with an emission wavelength of 450 nm and an excitation wavelength of 360 nm.

Molecular analysis

Seawater (2 L) collected from each sampling point in the mesocosm experiment was filtered through 0.2 µM pore cellulose acetate membrane filter (Sartorius Stedim Biotech, Göttingen, Germany). Nucleic acid was subsequently extracted from each of the filtered samples using the PowerWater DNA Isolation Kit (MoBio, Carlsbad, USA), according to the manufacturer's instructions. DNA was sent to MR DNA (www.mrdnalab.com, Shallowater, USA) for paired-end sequencing (2 x 300bp) of the 16S rRNA gene and the internal transcribed spacer (ITS) 1 rDNA variable region on the Illumina MiSeq platform (Illumina, San Diego, CA, USA), as detailed previously [START_REF] Phoma | Agulhas Current properties shape microbial community diversity and potential functionality[END_REF][START_REF] Vargas-Gastélum | Targeted ITS1 sequencing unravels the mycodiversity of deep-sea sediments from the Gulf of Mexico[END_REF]. The marine specific primer pair 515F-Y and 926R [START_REF] Walters | Improved bacterial 16S rRNA gene (V4 and V4-5) and fungal internal transcribed spacer marker gene primers for microbial community surveys[END_REF] was used to PCR amplify hypervariable regions V4 and V5 of the bacterial and archaeal 16S rRNA gene.

ITS region 1 was targeted and PCR amplified using the primer pair ITS1F [START_REF] Gardes | ITS primers with enhanced specificity for basidiomycetesapplication to the identification of mycorrhizae and rusts[END_REF] and ITS2 [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF]. All samples were barcoded on the forward primer to allow for multiplex sequencing.

Sequences were analysed using QIIME v 1.8.0 [START_REF] Caporaso | QIIME allows analysis of high-throughput community sequencing data[END_REF]. Sequences were demultiplexed according to barcode identity and were quality filtered according to QIIME default parameters, with exceptions for the ITS dataset. For the ITS dataset, the p score applied was altered to 0.45 to ensure that all reads >200bp were included in the dataset, provided that they met the other default quality parameters (i.e. q 19, r 3, etc.). Chimeras were detected both de novo and against the Greengenes reference library [START_REF] Desantis | Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB[END_REF] for 16S rRNA genes and the UNITE reference database for the ITS region [START_REF] Kõljalg | Towards a unified paradigm for sequence-based identification of fungi[END_REF] using usearch61 [START_REF] Edgar | UCHIME improves sensitivity and speed of chimera detection[END_REF]. Operational taxonomic units (OTUs) were clustered at a similarity threshold of 97% using the uclust method and the Greengenes and UNITE reference database for 16S rRNA gene and ITS region datasets, respectively. 16S rRNA gene OTUs were assigned a taxonomy using the Greengenes [START_REF] Desantis | Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB[END_REF] reference database. Taxonomy was assigned to ITS region OTUs using a combination of the UNITE database in QIIME and the Warcup Fungal ITS training set in the RDP classifier [START_REF] Wang | Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy[END_REF][START_REF] Kõljalg | Towards a unified paradigm for sequence-based identification of fungi[END_REF][START_REF] Deshpande | Fungal identification using a Bayesian classifier and the Warcup training set of internal transcribed spacer sequences[END_REF]. For the RDP classifier, taxonomic classifications were kept at the default confidence threshold of equal to or greater than 80% for classification to the root only, after which the default confidence threshold was decreased to 50% for all subsequent assignments. Following taxonomy assignment, singletons were removed from the datasets and both datasets were rarefied to the lowest number of sequences per sample (53,267 for 16S rRNA gene and 10,370 for ITS region), using a single rarefaction. All sequences are available on the SRA database under the following accession number: SRP126760

Network analysis

Ecological association networks were constructed using Random Matrix Theory (RMT) through the Molecular Ecological Network Analysis (MENA) pipeline (http://ieg2.ou.edu/MENA/) [START_REF] Deng | Molecular ecological network analyses[END_REF]. Individual networks were constructed for the three treatment groups (i.e. one low, medium, control) where all days were combined, excluding day zero (i.e. nine biological replicates per network, see Supporting Information for details). All networks were visualised in Cytoscape v 3.4.0 [START_REF] Shannon | Cytoscape: a software environment for integrated models of biomolecular interaction networks[END_REF].

Statistical analysis

All statistical analyses were performed using R statistical software version 3.2.2. (R Core Team, 2013). Alpha diversity was calculated using the vegan package [START_REF] Oksanen | Vegan: community ecology package version 2[END_REF]. Statistically significant differences in diversity measures between each treatment group were determined via the aov function (analysis of variance, ANOVA) in vegan and the post-hoc TukeyHSD test. Microbial community structure was visualised by non-metric multidimensional scaling (nMDS) ordination plots with the metamds function of vegan using Bray-Curtis dissimilarity and visualised using ggplot2. Significant differences in microbial community structure were determined via PERMANOVA (permutational analysis of variance) using the adonis and betadisper functions in vegan. The aov and tukeyHSD functions of vegan were used to determine significant differences between treatment groups for each of the enzyme activities. Correlations between taxon abundances, enzyme activities and pH were calculated using Spearman's rho, via the rcorr function of the Hmisc package.

Relationships between pH and enzyme activities, bacterial, archaeal and fungal diversity and taxon abundance were modelled using linear regressions (lm function). Fungal diversity relationships with pH were non-significant (p-value > 0.5), and therefore no regression line is shown. Topological roles of OTUs are indicated either as peripherals (Zi ≤ 2.5, Pi ≤ 0.62), connectors (Zi ≤ 2.5, Pi > 0.62), module hubs (Zi > 2.5, Pi ≤ 0.62) and network hubs (Zi > 2.5, Pi > 0.62). The abundance of peripherals, connectors and module hubs between the three treatment groups is indicated in (b), where the control group is the centremost donut ring, the medium group the middle ring and the low group the outer donut ring.

Figures legends

Fig. 1 :

 1 Fig. 1: Relationships between bacterial and fungal diversity and experimentally induced

Fig. 2 :

 2 Fig. 2: Non-metric multi-dimensional scaling (nMDS) ordination plots of bacterioplankton

Fig. 3 :

 3 Fig. 3: Relationship between extracellular enzymatic activities and experimentally induced

Fig. 4 :

 4 Fig. 4: Ecological bacterioplankton association networks for the control, medium and low pH

Fig. 5 :

 5 Fig. 5: Plot of the within-module connectivity (Zi) and between-module connectivity (Pi) of

Fig. 6 :

 6 Fig. 6: Conceptual diagram depicting ocean acidification effects on SO surface microbial

ACKNOWLEDGEMENTS

We are grateful to the National Research Foundation (NRF) (Grant ID 100052 SZdS) the South African National Antarctic Programme (SANAP 110717), and the University of Pretoria for funding. TPM also wishes to acknowledge the Fulbright Visiting Scholar Program for providing sabbatical funding. We wish to acknowledge Mr. A. van der Walt and Mr. M. Clayton for useful discussions on the content of the manuscript.


Johanna Zoppi

Jean-François Guillaume

Michel Neunlist

Samuel Chaffron
email: samuel.chaffron@univ-nantes.fr

MiBiOmics: An interactive web application for multi-omics data exploration and integration

Keywords: Multi-omics, Ordination, Biological networks, Data integration, R Shiny

Background: Multi-omics experimental approaches are becoming common practice in biological and medical sciences underlying the need to design new integrative techniques and applications to enable the holistic characterization of biological systems. The integrative analysis of heterogeneous datasets generally allows us to acquire additional insights and generate novel hypotheses about a given biological system. However, it can often become challenging given the large size of omics datasets and the diversity of existing techniques. Moreover, visualization tools for interpretation are usually non-accessible to biologists without programming skills.

Results: Here, we present MiBiOmics, a web-based and standalone application that facilitates multi-omics data visualization, exploration, integration, and analysis by providing easy access to dedicated and interactive protocols. It implements advanced ordination techniques and the inference of omics-based (multi-layer) networks to mine complex biological systems, and identify robust biomarkers linked to specific contextual parameters or biological states. Conclusions: Through an intuitive and interactive interface, MiBiOmics provides easy-access to ordination techniques and to a network-based approach for integrative multi-omics analyses. MiBiOmics is currently available as a Shiny app at https://shiny-bird.

Background

The holistic characterization of biological systems is extending our knowledge about the functioning of organisms and natural ecosystems. Today, their multi-omics characterization is becoming standard, thus novel methodologies and easily accessible tools are required to facilitate the study of associations and interactions within and across omics layers (e.g. (meta-)genome, (meta-)transcriptome, metabolome) and scales (e.g. cells, organs, holobionts, communities). The analysis of single omics datasets has helped to identify molecular signatures associated to phenotypes of interest. However, it usually does not allow to predict mechanisms underlying phenotypic variabilities. Although multi-omics information is not sufficient to identify causes and consequences of a biological process, it can contribute to delineate key players sustaining it [START_REF] Hasin | Multi-omics approaches to disease[END_REF]. Indeed, exploring a biological system across several omics layers enable to capture additional sources of variability associated with a variation of interest and potentially to infer the sequence of events leading to a specific process or state [START_REF] Li | A review on machine learning principles for multi-view biological data integration[END_REF]. Within the last decade, multi-omics integrative approaches have been applied across various fields including microbial ecology [START_REF] Heintz-Buschart | Integrated multi-omics of the human gut microbiome in a case study of familial type 1 diabetes[END_REF], genetics [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in late-onset Alzheimer's disease[END_REF] and personalized medicine [START_REF] Chen | Personal omics profiling reveals dynamic molecular and medical phenotypes[END_REF]. As of today, several integrative methods have been developed, but are often specific to a given experimental design, data type or a precise biological question [START_REF] Paliy | Application of multivariate statistical techniques in microbial ecology[END_REF]. Considering the multiplicity of existing techniques, the selection of an appropriate workflow is challenging for biologists, especially when it comes to the representation of several system-level omics layers and its interpretation. There is a clear need for accessible (web) tools to facilitate the integration, analysis and representation of multi-omics datasets through an intuitive and guided approach. MiBiOmics aims to provide established and novel techniques to reveal robust signatures in high dimensional datasets [START_REF] Guidi | Plankton networks driving carbon export in the oligotrophic[END_REF] through a graphical user interface allowing to perform widely applicable multi-omics analyses for the detection and description of associations across omics layers. Available as a web-based and a stand-alone application, it gives access to several R packages and tools to help users who are not familiar with programming to load and explore their data in a simple and intuitive way. MiBiOmics allows the parallel study of up to three omics datasets, as well as the in-depth exploration of each single dataset. It also provides easy access to exploratory ordination techniques and to the inference of (multi-layer) correlation networks enabling useful dimensionality reduction and association to contextual parameters. The user can then compare results from these different approaches and cross-validate multi-omics signatures to generate confident novel hypotheses.

Implementation

MiBiOmics is implemented in R (Version 3.6.0) as a Shiny app providing an interactive interface to perform each step of a single-or multi-omics data analysis (Figure 1). MiBiOmics is also accessible as a standalone application that can be easily installed via Conda (Version 4.6.12). The application is divided into five sections as described below:

Data Upload

Within MiBiOmics, the user can upload up to three omics datasets, allowing the data exploration and network analysis of a single-or multi-omics dataset. There must be common samples between omics datasets in order to perform all analyses provided by the application. An annotation table describing external parameters (e.g. pH, site of extraction, physiological measures) needs to be provided. These parameters may be quantitative or qualitative, and available for each sample. An additional taxonomic annotations table can be uploaded when one omics table corresponds to microbial lineages (e.g. as Operational Taxonomic Units (OTUs) or Amplicon Sampling Variants (ASVs)). Following data upload, the user can filter, normalize and transform each data matrix using common methods, such as the center log ratio (CLR) transformation to deal with the compositional nature of sequencing data, or filtration based on prevalence. In this section, it is also possible to detect and remove potential outlier samples. To allow new users to easily test the functionality of MiBiOmics, we provide two example datasets: the breast TGCA datasets from The Cancer Genome Atlas [8] allows to explore assocations between miRNAs, mRNAs and proteins in different breast cancer subtypes; and a dataset from the Tara Oceans Expeditions [START_REF] Sunagawa | Ocean plankton. Structure and function of the global ocean microbiome[END_REF][START_REF] Mariette | Unsupervised multiple kernel learning for heterogeneous data integration[END_REF] to explore prokaryotic community compositions across depth and geographic locations.

Data Exploration

In this section, several ordination plots (Principal Component Analysis (PCA), Principal COordinates Analysis (PCoA)) [START_REF] Dixon | VEGAN, a package of R functions for community ecology[END_REF] are dynamically produced to visualize and explore relationships between samples, and to identify main axes of variation in each dataset. When OTUs or ASVs are uploaded with their taxonomic annotations, it is possible to obtain a relative abundance plot describing the proportion of lineages at a given taxonomic level (e.g. Phylum, Family, Genus or Species) in each sample.

Network Inference

The network inference section allows to perform a Weighted Gene Correlation Network Analysis (WGCNA [START_REF] Langfelder | WGCNA: An R package for weighted correlation network analysis[END_REF]). Help sections are available to assist the user with parametrization, notably for optimizing the scale-free topology of the network. Here, WGCNA networks can be inferred for each uploaded omics dataset. We strongly advise users to read the WGCNA original publication and associated tutorials for this step of the analysis.

Network Exploration

The network exploration section allows to compute and explore significant associations between subnetworks or modules (e.g. of genes, transcripts, metabolites), and communities (of lineages) delineated from each omics layer, which are denser and contain highly correlated features. Each module is associated to all external parameters provided in the annotation table and correlations are visualized as a heatmap (Figure 2A). Modules associated to parameters of interest can be further analyzed. The user can also identify which samples are contributing the most to the delineation of a specific module (Figure 2B). In case an OTUs/ASVs table is provided with taxonomic annotations, the relative abundance of lineages contributing to each module can be visualized as barplots. In addition, OPLS (Orthogonal Partial Least Square) regressions [START_REF] Wehrens | The pls Package: Principal Component and Partial Least Squares Regression in R[END_REF] can be performed using a selected module components as features in order to estimate its capacity to predict a given contextual parameter, and are useful to cross-validate a module-parameter association. The results of this analysis are represented as hive plots with two axes. On the x-axis, the module features are ordered according to their Variable Importance Projection (VIP) score (a measure of their weight in the OPLS regression), while on the y-axis they are ordered according to their correlations to an external parameter of interest (Figure 2C).

Multi-Omics Analysis

Here, MiBiOmics allows users to detect and study associations across omics datasets. Ordination methods including Procrustes analysis [START_REF] Dixon | VEGAN, a package of R functions for community ecology[END_REF] and multiple coinertia [START_REF] Meng | A multivariate approach to the integration of multi-omics datasets[END_REF] are useful to compute and visualize the main axes of covariance, and also to extract multi-omics features driving this covariance. This central section of MiBiOmics implements an innovative approach for detecting robust links between omics layers. All modules delineated within each omics-specific network are associated to each other by directly correlating their eigenvectors. Here, the dimensionality reduction of each omics dataset through module definition ensures a small number of correlations, thereby increasing the statistical power for detecting significant associations between omics layers. For visualization, a hive plot helps summarizing significant associations between each module as a multi-layer network integrating links between omics-specific modules as well as their association to contextual parameters (traits or phenotypic characteristics). In this hive plot, each axis represents the network of a given omics layer. Corresponding modules are ordered on the axes according to their association to a contextual parameter of interest selected by the user. Modules with no significant associations are not depicted. Significant associations between omics-specific modules are represented, and individual associations between modules can also be visualized as heatmaps and dataframe. Conveniently, the user can also select modules of interest to investigate pairwise correlations between modules' features and delineates group of modules associated together and to an external parameter. Importantly, all figures generated by the application (PCA, PCoA, relative abundance plots, WGCNA outputs, hive plots, multiple co-inertia, Procrustes plots, correlograms, bipartite networks) can be downloaded (as svg or pdf files), and if aWpplicable, supplementary information provided as a csv file (WGCNA modules information, eigenvalues and coinertia drivers).

Results and Discussion

MiBiOmics enables the exploration, integration, analysis and visualization of up to three omics datasets. Through the primary exploration of a dataset, the inference of biological networks and the extraction of multi-omics associated features, the application provides a ready-to-use analysis pipeline to interactively explore sources of variability and variables of interest in a given biological dataset, as well as associations between multi-omics features in holistic studies. The inference of networks from omics features is useful to represent and model the complex architecture of putative interactions in biological systems. In addition, networks provide a way to reduce the dimensionality of a dataset by delineating cohesive groups of co-varying, often functionally related features, that can then be associated to contextual or phenotypic characteristics of interest [START_REF] Hasin | Multi-omics approaches to disease[END_REF]. A key functionality of MiBiOmics is the multi-omics adaptation of WGCNA [START_REF] Langfelder | WGCNA: An R package for weighted correlation network analysis[END_REF] to explore association across omics datasets via a network-based approach. As shown in Figure 2A, the interface provides the ability to interactively probe associations across omics layers of different breast cancer subtypes [8] within each network and their association to patient parameters. The original WGCNA outputs are provided by the application to deepen the analysis between modules and external parameters (Figure 2B). In addition, we provide the user with the possibility to perform an OPLS regression for modules of interest to evaluate the robustness of these variables to predict a given trait or phenotype. Figure 2C is an example of an OPLS regression using WGCNA module variables as features. On the x-axis the features of the red module are ordered according to their VIP score (their importance for the module), and on the y-axis according to their correlation to the subtype parameter. This figure highlights how central features of a WGCNA module relate to an external parameter. The exploratory multi-omics analysis allows to study the main axes of covariance across omics profiles and give the ability to discover and select variables implicated in an association between omics datasets. The concomitant application of (multiple) co-inertia (Figure 3A) and/or Procrustes ordination techniques with the exploration of multi-omics correlations between WGCNA modules of distinct omics layer (Figure 3B), provides a complementary vision of multi-omics relationships. The MiBiOmics interface allows to explore WGCNA modules of interest to directly infer significant associations between features from distinct omics layers (Figure 3C). In a multi-omics adaptation, WGCNA can be used to delineate a group of modules associated together and to a parameter to interest and extract features of different omics nature but related to each other. While an interactive version of WGCNA already exists [START_REF] Sundararajan | Shiny-Seq: advanced guided transcriptome analysis[END_REF], MiBiOmics goes beyond by providing a multi-omics strategy to identify correlated modules across omics layers and generate novel hypotheses. Associating modules across different datasets has already been performed in the original WGCNA article [START_REF] Langfelder | WGCNA: An R package for weighted correlation network analysis[END_REF] and reproduced in several studies. For example, the overlap of modules between transcriptional profiles of different tissue [START_REF] Xiao | Multi-tissue Analysis of Co-expression Networks by Higher-Order Generalized Singular Value Decomposition Identifies Functionally Coherent Transcriptional Modules[END_REF] was assessed, as well as a comparison between proteomics and gene expression profile of modules in a cohort of alzheimer patients [START_REF] Seyfried | A Multi-network Approach Identifies Protein-Specific Co-expression in Asymptomatic and Symptomatic Alzheimer's Disease[END_REF]. In both cases, the association between modules was determined by overlapping identical features (e.g. same genes in a given reference genome) within each module, a method which is not applicable when omics datasets do not contain similar data types or refer to the same biological system. In MiBiOmics, we enable the inference of relationships between omics layers within an entire biological system (e.g. holobiont) or ecosystem (e.g. the plankton), which makes it more widely applicable and especially suited for omics-based environmental studies. While there is a growing need for multi-omics analysis tools, a few applications have been developed to fulfill the need of integrating and mining these complex and heterogeneous datasets. Several methods such as MONGKIE [START_REF] Jang | MONGKIE: An integrated tool for network analysis and visualization for multi-omics data[END_REF], are based on prior knowledge and integrate data by projecting them on known metabolic networks and biological pathways. More generally, existing multi-omics pipelines are focusing on certain data types (Metabolomics with MetaboAnalyst [START_REF] Chong | MetaboAnalyst 4.0: Towards more transparent and integrative metabolomics analysis[END_REF]) or on disease-related mechanisms (MergeOmics [START_REF] Shu | Mergeomics: Multidimensional data integration to identify pathogenic perturbations to biological systems[END_REF]). More widely applicable methods exist, such as the R package mixOmics [START_REF] Rohart | mixOmics: An R package for 'omics feature selection and multiple data integration[END_REF] that provides several semi-supervised methodologies often based on ordination techniques. We compared methods integrated in MiBiOmics (see Supplementary information for details) to the mixOmics DIABLO methodology [START_REF] Singh | DIABLO: an integrative approach for identifying key molecular drivers from multi-omics assays[END_REF]. Only few multi-omics features associated to breast cancer subtype in the TGCA dataset were extracted by all three methods (n=32, Figure 4A and Supplementary table 1). Both methods integrated in MiBiOmics (i.e. multiple co-inertia and multi-omics WGCNA) and DIABLO extracted mostly distinct features (Figure 4A) underlying the probable complementarity of these multi-omics integrative strategies. Scores attributed by each method to the common set of extracted features were also dissimilar (Figure 4BCD and Supplementary table 1). This may be explained by the fact that these methods implement fundamentally different approaches to features extraction and selection, which confirms the complementary nature of each analysis. For comparing the predictive power of models integrating features extracted by each method, we performed Sparse Partial Least Square Discriminant Analysis (sPLS-DA) and computed the corresponding mean AUC scores (Figure 4A and Supplementary Figure 1). All models can be considered to be highly predictive of the cancer subtype phenotypes with the miBiOmics multiomics WGCNA methodology obtaining the highest AUC score (AUC=0.9952), while the multiple coinertia analysis performed very well too (AUC=0.9903). Features extracted by the DIABLO method from mixOmics obtained the lowest score (AUC=0.9728) but remained highly predictive. Generally, these methods may benefit from an enrichment method applied to the list of extracted drivers [START_REF] Meng | A multivariate approach to the integration of multi-omics datasets[END_REF]. Overall, MiBiOmics provides two complementary methods to extract associated variables between omics layers and in relationship with a trait of interest. While it is well suited to generate new hypotheses about molecular processes, it can not infer causal mechanisms between omics features and phenotypes, which would require experimental validations but that can actually be guided by MiBiOmics results. In addition to provide an exploratory and integrative framework for multi-omics studies, MiBiOmics distinguishes itself by providing a powerful dimensionality reduction and unsupervized method combining both ordination and graph-based techniques, which enables to study complex biological systems as a whole. Importantly, it also integrates contextual information by linking multi-omics signatures to qualitative and quantitave contextual parameters.

Conclusion

MiBiOmics is an interactive web-based (and standalone) application to easily and dynamically explore associations across omics datasets. Through an innovative network-based integrating strategy, it can help biologists to identify potential mechanisms of interactions and generate novel hypotheses. The core of the application lies behind the reduction of dimensionality across omics datasets to efficiently link them at the molecular level, and to identify biomarkers associated with a given trait or phenotype. The MiBiOmics pipeline facilitates the exploration, integration and analysis of multi-omics datasets to a broad audience by providing scientists a powerful way to predict and explore putative molecular mechanisms underlying complex phenotypes across a wide range of biological scales and systems.

Availability and requirements

Project name: MiBiOmics Project home page: https://gitlab.univ-nantes.fr/combi-ls2n/mibiomics Operating system(s): Platform independent Programming language: R Other requirements: for the local installation Conda 4.6.12 or Docker. License: AGPL-3 Any restrictions to use by non-academics: No restrictions.

Figures

Figure 1 MiBiOmics framework. The MiBiOmics workflow can be divided into three main tasks: Data preprocessing, data exploration, and multi-omics integration. The data preprocessing task is dedicated to data upload, data filtration, normalization and transformation. The data exploration task implements classical clustering methods, PCA, PCoA and WGCNA correlation networks that can be applied to each omics dataset separately. Finally, the multi-omics integration task allows the user to perform multi-omics exploration, integration and analyses using ordination techniques (multiple co-inertia and Procrustes analysis), and multi-omics network inference.

Figure 2

 2 Figure 2 MiBiOmics networks exploration. MiBiOmics networks analysis visualizations (network exploration section) from the analysis of the The Cancer Genome Atlas Network [8] breast TGCA datasets. A. Correlation heatmap displaying associations of interest between WGCNA modules and contextual parameters. B. The upper panel indicates the contribution of each sample in the red module delineation. The lower panel indicates the corresponding subtype value for each sample. Here, Basal samples positively contributes module red, while Her2, LumA and LumB negatively contributes to module red. C. Hive plot displaying the red module's features according to their VIP scores, correlations to the subtype parameter and their relationships.

Figure 3

 3 Figure 3 MiBiOmics multi-omics integration. MiBiOmics multi-omics integration visualizations (multi-omics integration section) from the analysis of the The Cancer Genome Atlas Network [8] breast TGCA datasets. A. A multiple co-inertia plot integrating 3 omics layers and extracted miRNA, mRNA and protein drivers. B. Hive plot displaying modules of each omics network and their associations. C. Bipartite network between mRNA features of the red mRNA module and miRNA features of the red miRNA module (Spearman Correlation > 0.35).

Figure 4

 4 Figure 4 MiBiOmics and mixOmics comparison. A. Venn diagram displaying common and distinct features extracted by DIABLO (mixOmics), multiple coinertia and multi-omics WGCNA (MiBiOmics). Area Under the Curve (AUC) scores were computed to compare features selection and model performance of each method. B. Selected features' weights comparison between multiple coinertia and multi-WGCNA on their common subset of features. C. Selected features' weights comparison between multiple coinertia and DIABLO (mixOmics) on their common subset of features. D. Selected features' weights comparison between DIABLO (mixOmics) and multi-WGCNA method on their common subset of features.

Acknowledgements

The authors wish to thank Mélanie Fouesnard, Damien Eveillard, Philippe Bordron, Gaëlle Boudry, Catherine Michel and Simon Beck for their useful feedback while testing and using the application. We also thank the bioinformatics core facility of Nantes (BiRD -Biogenouest) for providing computing resources and support.

Funding

This work has received financial support from the Region Pays de la Loire (MiBioGate 2016-11179 to M.N.), the French National Institute of Health and Medical Research and the CNRS through the MITI interdisciplinary program Modélisation du Vivant [GOBITMAP to S.C.].

univ-nantes.fr/app/Mibiomics and as a

Availability of data and materials

The datasets provided as example within MiBiOmics application are available in the data repository, at https://gitlab.univ-nantes.fr/combi-ls2n/mibiomics .

Competing interests

The authors declare that they have no competing interests.

Authors' contribution:

JZ and SC participated equally in the design and development of MiBiOmics as a standalone web-application, and writing the manuscript. MN contributed to the writing of the manuscript. JFG was responsible for the distribution and maintenance of MiBiOmics web-server at https://shiny-bird.univ-nantes.fr/app/Mibiomics. All authors read and approved the final version of the manuscript.

Author details 1 Université de Nantes, INSERM, TENS, Nantes, France. 2 Université de Nantes, CHU Nantes, Inserm, CNRS, SFR Santé, Inserm UMS016, CNRS UMS 3556, 44000 Nantes, France. 3 Université de Nantes, CNRS UMR 6004, LS2N, 44000 Nantes, France. 4 Research Federation (FR2022) Tara Oceans GO-SEE, Paris, France.