Bacterioplankton community resilience to ocean acidification : evidence from microbial network analysis

Yu Wang1‡, Rui Zhang1‡, Qiang Zheng1, Ye Deng2, Joy D. Van Nostrand3, Jizhong Zhou3,4,5*, and Nianzhi Jiao1* State Key Laboratory of Marine Environmental Science, Institute of Marine Microbes and Ecospheres, Xiamen University, Xiamen, Fujian 361005, People’s Republic of China CAS Key Laboratory of Environmental Biotechnology, Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing 100085, People’s Republic of China Institute for Environmental Genomics and Department of Microbiology and Plant Biology, University of Oklahoma, Norman, OK 73019, USA State Key Joint Laboratory of Environment Simulation and Pollution Control, School of Environment, Tsinghua University, Beijing 100084, China Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA


Introduction
Since the industrial revolution, the impact of human activity on the global climate has increased greatly as a result of increasing carbon dioxide (CO 2 ) emissions from anthropogenic sources.The uptake of anthropogenic carbon dioxide by the ocean has caused a decrease of pH 0.1 units (ocean acidification, OA; IPCC, 2015).Previous studies have demonstrated that some phototrophic communities, like sea grass (Zimmerman et al., 1997;Jiang et al., 2010), diatoms (Riebesell et al., 1993;Baragi et al., 2015;Taucher et al., 2015), and coccolithophorids (Hiwatari et al., 1995;Lv et al., 2015), have higher photosynthesis rates under a higher partial pressure of carbon dioxide (pCO 2 ).However, the response of marine bacterioplankton, a crucial player in marine biogeochemical cycling (Azam, 1998;Jiao et al., 2010), to OA is not well understood at present.Joint et al. (2011) proposed that "marine microbes possess the flexibility to accommodate pH change".Indeed, several mesocosm studies have found that elevated pCO 2 has little effect on bacterial communities in the Arctic Ocean (Allgaier et al., 2008;Tanaka et al., 2008;Oliver et al., 2014).The abundance and activity of bacteria in these communities did not differ statistically under three different pCO 2 treatments (350, 700, and 1085 matm;Allgaier et al., 2008).Furthermore, phylogenetic diversity analysis revealed no clear effect of elevated pCO 2 on a bacterioplankton assemblage in the high Arctic Ocean (Monier et al., 2014).In contrast, a few studies have demonstrated that elevated pCO 2 has some influence on microbial community composition (Krause et al., 2012;Bowen et al., 2013).Other studies have also shown that the production and growth rates of bacterial isolates were strongly affected by high pCO 2 and low pH (Takeuchi et al., 1997;Labare et al., 2010;Teira et al., 2012).For example, Vibrio alginolyticus growth was suppressed under low pCO 2 levels (Michael et al., 2010).In contrast, high pCO 2 stimulated the growth efficiency of one Flavobacteria strain (Teira et al., 2012).Moreover, the rate of microbial ammonia oxidation is inhibited by reduced pH in both surface and deep seawater (Huesemann et al., 2002).Conflicting results from the population/ecosystem and species levels indicate that community may play a crucial role in determining the response of microbes to OA.In addition to community structure and the number of species, micro-organism interaction is an important component of diversity (Olesen et al., 2007).For example, ecological networks among different species bridge ecosystem complexity and stability (Montoya et al., 2006).The interaction between plants and pollination enhances the relative resistance of plants to environmental disturbance (Sole and Montoya, 2001;Memmott et al., 2004).Compared with network investigations among macro-organisms (Elton, 1927;Paine, 1980;Bascompte et al., 2003), microbial interactions and ecological networks are understudied and have only been investigated recently with advances in molecular technology (Sherr and Sherr, 2008;Chaffron et al., 2010;Steele et al., 2011).Several microbial network investigations have been conducted on temporal variation in microbial communities and the interactions among bacteria, phage, and protists at the San Pedro Ocean Time-series (SPOT) station off southern California (Chow et al., 2013(Chow et al., , 2014)).Bacterial interaction has also been proposed as a determinant of phytoplankton bloom dynamics (Tan et al., 2015).Currently, only a few investigations have been performed to illustrate the interaction of microbial communities in response to climate change in soil environments (Zhou et al., 2010;Tu et al., 2015), but no studies have addressed interactions among microorganisms under elevated pCO 2 in the marine environment.
A large-scale mesocosm experiment, European Project on OCean Acidification (EPOCA), was set up in the pelagic coastal area of the Arctic Ocean to assess the impact of OA on the marine ecosystem.Previous EPOCA studies have demonstrated that the phytoplankton community was more susceptible to elevated pCO 2 than the bacterioplankton community, where pico-phytoplankton and diatoms were stimulated by nutrient addition and elevated pCO 2 , resulting in increased primary production (Riebesell et al., 2013b).Bacterial activity was subsequently spurred by stimulated primary production and elevated pCO 2 , leading to a decrease in net production (Engel et al., 2013;Motegi et al., 2013).However, bacterial production and respiration were generally similar among the pCO 2 mesocosms during the incubation period (Motegi et al., 2013).Nonetheless, bacterioplankton activity was closely related to phytoplankton under both ambient and elevated pCO 2 , and limited by both nutrient availability and elevated pCO 2 (Riebesell et al., 2013b).The general diversity and composition of microbial communities has been evaluated by Automated Ribosomal Spacer Analysis (Sperling et al., 2013), T-RFLP, clone libraries (Zhang et al., 2013), and next-generation sequencing (Roy et al., 2013), all of which indicated that increased pCO 2 did not significantly impact the communities.However, because of limited sample numbers and sequencing depth, the above dataset cannot be used to explore the interactions among bacterial groups under elevated pCO 2 .Thus, we examined the structure and diversity of microbial communities at this site using Illumina sequencing and performed phylogenetic molecular ecological network (pMEN) analysis using the Molecular Ecological Network Analysis Pipeline (MENAP; Deng et al., 2012) for samples collected from nine mesocosms, each having a different pCO 2 level (175-1085 matm), which simulated different acidification scenarios.We tested the null hypothesis that elevated pCO 2 would have an insignificant or indirect impact on bacterial assemblages as proposed by Joint et al. (2011), with a focus on the interactions among bacterioplankton under different pCO 2 concentrations.Our results indicated an insignificant effect of elevated pCO 2 on microbial assemblage interactions and connectivity and provide a deep insight into our understanding of bacterioplankton assemblage resistance to climate change.

Experimental set-up and sampling
A large-scale mesocosm experiment, supported by EPOCA, was conducted in the Arctic Ocean at Kongsfjorden, Svalbard, Norway (78856.2′ N, 11854.6 ′ E), from June to July 2010.Experimental set-up details have been described previously (Riebesell et al., 2013a).Briefly, nine mesocosms (45 m 3 volume), each simulating a different pCO 2 concentration ( 175,180,250,340,450,600,675,860, and 1085 matm) by bubbled CO 2 , were deployed into the Arctic Ocean.On day 13, nutrients were added to induce phytoplankton bloom development (Schulz et al., 2013).
Sample collection was carried out according to Zhang et al. (2013).Briefly, seawater (2 l) collected from nine mesocosms was filtered onto a 0.22 mm pore filter (GTTP, Millipore) to collect microbial cells.The filters were stored at 2808C until DNA extraction.In total, we obtained 165 environmental DNA samples from the nine mesocosms over 19 sampling points (not all treatments contain all time points; Supplementary Table S1).Geochemistry and biological parameters were downloaded from the EPOCA website (http:// www.epoca-project.eu;Nisumaa et al., 2010).The geochemistry parameters included temperature, salinity, pH, oxygen, particle organic carbon (POC), particle organic nitrogen, particle organic phosphorus, NO 3 , Si, PO 4 , and NH 4 ; the biological parameters included bacterial production, chlorophyll a concentration, bacterial abundance, bacterial respiration, bacterial biomass production, viral abundance, and pico-phytoplankton abundance.
DNA extraction, Illumina high-throughput sequencing, and analysis DNA was extracted by a freeze-grinding method as described previously (Zhang et al., 2013).The forward (515F, 5 ′ -GTGCCAGCM GCCGCGG-3 ′ ) and reverse (806R, 5 ′ -GGACTACHVGGGTWTC TAAT-3 ′ ) primers were used to target the bacterial 16S rDNA gene V4 hyper-variable region (Caporaso et al., 2012).PCR amplification was performed as described by Caporaso et al. (2010).The amplicons were then paired-end sequenced using the Illumina MiSeq platform (Illumina, Inc., San Diego, CA, USA) and the reads were then analysed through an in-house sequence analysis pipeline (IEG sequence analysis pipeline, http://zhoulab5.rccc.ou.edu).Briefly, the sequence length was trimmed (minimum length, 100 bp), Btrim was used to remove sequencing adaptors and low quality regions, and sequences with quality scores ,20 were then removed 866 Y. Wang et al. (Kong, 2011).Thereafter, UCHIME was employed to remove the chimeric sequences (Edgar et al., 2011).The number of sequences from each sample was resampled randomly to 35 000.The sequences were then clustered using the CD-HIT algorithm (Li and Godzik, 2006).The taxonomic assignment was determined using the RDP (Ribosome Database Project) classifier against the RDP 16S rRNA database based on BLAST (Wang et al., 2007).

pMEN construction and network analysis
Sequence alignment was carried out by the PyNASTalignment method against pre-aligned reference 16S rRNA sequences (Caporaso et al., 2009), followed by approximately maximum-likelihood tree construction using the FastTree method with default settings (Price et al., 2009;Morgan et al., 2010).Comparison of microbial community diversity was performed using the Unifrac method in a phylogenetic context (Lozupone andKnight, 2005, 2007;Lozupone et al., 2006) followed by principal coordinates analysis (PCoA).pMENs analysis was performed based on the relative abundance of all samples through the molecular ecological networks analysis pipeline (http://ieg2.ou.edu/MENA; Zhou et al., 2010Zhou et al., ,2011;;Deng et al., 2012).The whole process and details are given in a previous MENA study (Deng et al., 2012); we briefly introduced the major process in the present study.OTUs that appeared in 13 or fewer samples were removed.The relative abundance of OTUs was log-transformed and missing values were filled with a very small number (0.01) if paired valid values were available.This step ensured a more statistically reliable correlation coefficient between two OTUs.Following data preparation, the relative abundance of each OTU in each sample was used to generate the similarity matrix as the foundation for subsequent steps.Similarity matrices (adjacency matrix) were created for each network based on the pairwise Pearson correlation coefficient across the time-series (one time points lag).To obtain a confident network construction, the threshold of pairwise Pearson correlation coefficient values between OTUs was identified by an RMT-based approach that observed a transition point of nearest-neighbour spacing distribution of eigenvalues from Gaussian to Poisson distribution (Zhou et al., 2010).To compare networks, an identification cut-off of 0.84 (0.85 for 175 and 860 matm treatments) was used to construct the microbial community networks.Additionally, ecological networks predicted by R 2 (R 2 .0.8) generated based on the random matrix theory (RMT) should be scale-free (Zhou et al., 2010).An identification cut-off of 0.84 for 175 and 860 matm treatments did not obtain an R 2 .0.8, where 0.85 could.To evaluate whether or not the constructed networks were random, we employed a permutation-based null model analysis with 100 permutations and kept the number of nodes and links constant.A Student's t-test was performed to evaluate significance differences between random and experimental networks.Once the pMEN was determined, the topological indices were calculated based on the adjacency matrix.Module detection of each network was based on fast greedy modularity optimization (Newman, 2006a).Identification of key module members was based on within-module connectivity (Z i ) and among-module connectivity (P i ) of each node (Olesen et al., 2007).Eigen-gene, which is a linear combination of genes and eigenvalues (Deng et al., 2012), was detected by singular value decomposition analysis (Alter et al., 2000).Thereafter, trait-based gene significance (GS), which is the correlation of relative abundance of each OTU to a sample trait (e.g.pCO 2 , temperature), was calculated.Networks were visualized with the NetworkAnalyzer plugin of the Cytoscape 3.1.1(Assenov et al., 2008).Mantel tests were performed to identify relationships between pMENs and biological and chemical variables (Zhou et al., 2010).

Nucleotide sequence accession number
Sequences obtained in this study were uploaded to the MG-RAST database under the ID number 4626065.3.

Results
A previous study sequenced samples from this same experimental site, collected across nine time-points, using the Illumina HiSeq2000 platform (Roy et al., 2013); however, network analysis studies and our preliminary tests suggested that this sampling frequency does not result in accurate and reliable network analysis.Therefore, we expanded this sample set to include 19 time-points for each of nine pCO 2 treatments, nine more time-points than Roy et al. (2013;Supplementary Table S1).These DNA samples were also used to investigate the bacterioplankton community by T-RFLP in a previous study (Zhang et al., 2013).Results from different sequencing platforms or different runs on the same platform can affect the estimation of microbial diversity (Zhou et al., 2013), so all samples were sequenced together on the Illumina MiSeq platform.The microbial community structure revealed in our study (discussed below) did not show obvious differences from those of Roy et al. (2013).In total, 8 900 000 paired-end reads were generated after quality control in this study, with an average length of 250 bp (Supplementary Table S1).To increase the reliability of comparison among samples, the sequences were resampled with a threshold of 35 000.

Bacterial community composition and diversity
Similar bacterial community composition and structure were observed in the nine pCO 2 treatments during the entire incubation period (Table 1, Figure 1, and Supplementary Figure S2), which is consistent with previous results (Roy et al., 2013;Sperling et al., 2013;Zhang et al., 2013).The Proteobacteria, Bacteroidetes, Actinobacteria, and Cyanobacteria/chloroplast were predominant in all nine mesocosms.The relative abundance of these microorganisms varied over the incubation period.For example, the Proteobacteria dominated at the beginning of the incubation period and then decreased in relative abundance; while the Bacteroidetes exhibited a relatively high abundance at the end of the incubation, which could be a response to the phytoplankton bloom (Figure 1).Meanwhile, the rarefaction curves for all samples revealed variations in richness and diversity among samples (Supplementary Figure S1).The bacterioplankton community diversity, including Chao index and PD, as well as richness, decreased during the incubation period in all pCO 2 treatments (Supplementary Figures S2 and S3).The weighted PCoA results revealed that the bacterial community structures were altered by nutrient addition, although there was no clear pattern among different pCO 2 treatments (Supplementary Figure S4).

Microbial community interactions
The increased number of sequences obtained by high-throughput sequencing technology afforded us the unprecedented opportunity to investigate the interaction and connectivity within these microbial communities using network analysis.pMEN analysis is a novel RMT-based framework for studying microbial interaction, which has been used in long-term CO 2 experiments (Zhou et al., 2010(Zhou et al., , 2011;;Tu et al., 2015).In the pMEN analysis, a node in pMEN indicates an OTU, while a link indicates a correlation between two connected nodes.The connectivity, also called node degree, indicates the connection strength of a node; therefore, Bacterioplankton community resilience to OA higher average connectivity means a more complex network.Geodesic distance means the shortest path between two nodes, while a smaller average geodesic distance means the nodes in a network are closer.The average clustering coefficient means the extent of a module in a network, where the clustering coefficient describes how well a node is connected with its neighbours.The network can be separated into communities or modules (Newman, 2006b).The modularity value indicates how well a network can be divided into modules.In our study, elevated pCO 2 generally had no consistent influence on the connectivity of micro-organisms observed via pMENs analysis (Table 1).For example, the size (number of total nodes) of pMENs under 175 and 180 matm pCO 2 were smaller than for higher-level treatments (250, 425, 600, and 675 matm), while they were similar to those from the 860 and 1085 matm mesocosms.Additionally, the total number of links and average connectivity were similar among the nine pCO 2 treatments.Average geodesic distance and average clustering coefficients revealed a discrepancy change along the pCO 2 gradient.The modularity of nine pMENs, which contained at least four modules each, suggested that the bacterial communities in these mesocosms were highly complex (Olesen et al., 2007).Strong negative correlations were observed in all nine pMENs, while positive correlations primarily occurred in submodules (Figure 2).These results suggested potential competition among bacterial lineages in these ecosystems.
Connectivity analysis among or within the modules showed that different microbial clades played different roles in the pMENs (Figure 3).From an ecological perspective, the peripherals may represent specialists, whereas module hubs and connectors may be more generalists and network hubs may be super-generalists (Olesen et al., 2007;Deng et al., 2012).The relatively higher abundance of Alphaproteobacteria and Flavobacteria suggest that they play crucial roles in the pMENs.The number of connectors was greater under 250, 600, and 1085 matm pCO 2 , which suggested a greater complexity in these microbial communities (Supplementary Figure S5).Larger numbers of OTUs belonging to Alphaproteobacteria and Flavobacteria were observed under 250 and 600 pCO 2 .Network module hubs are another important factor in network topology.Compared with connectors, there were fewer module hubs, but they were dominated by Gammaproteobacteria (Supplementary Table S4).Unexpectedly, no network connectors or module hubs were detected under 175 and 860 matm.Consistent with previous studies, no network hub was detected in any pMEN (Zhou et al., 2010(Zhou et al., , 2011;;Deng et al., 2012).These results indicated that Alphaproteobacteria, Gammaproteobacteria, and Flavobacteria play an important role in maintaining the stability of microbial ecosystems in these mesocosms.
The Mantel test was performed to investigate the correlation between environmental variables and micro-organism interactions.The OTU significance (GS) in the pMENs were calculated, which were defined as the square of the Pearson correlation coefficient (r 2 ) of OTU abundance profile with environmental traits (Deng et al., 2012), and did not significantly correlate with either biological or chemical parameters in almost all of the pMENs, except for the 860 matm treatment (Supplementary Table S2).Furthermore, rare individual variables were significantly correlated with the network topological properties.However, the 250 and 340 matm pMENs were significantly correlated with particle organic phosphorus and NO 2 concentration, respectively.The 675 and 860 matm pMENs were significantly correlated with incubation day and pH, and bacterial abundance, viral abundance, and particle organic nitrogen, respectively.These significant correlations indicated connectivity  among these micro-organisms under 860 matm pCO 2 because both were significantly associated with changes in biological and geochemical variables.However, there was no clear trend in these correlations between topology properties and environmental variables along the pCO 2 gradient, which suggested that elevated pCO 2 had little effect on the global interaction properties.Mantel analysis was also used to investigate the correlation between environmental variables and pMEN microbial composition.Similarly, the Mantel test results between pMEN topology and environmental variables revealed that few OTUs were significantly correlated with environmental variables (Supplementary Table S3).Cyanobacteria/chloroplast were significantly correlated with environmental properties under the 180 and 1085 matm pCO 2 treatments.Betaproteobacteria and Gammaoproteobacteria were significant in the 675 matm pCO 2 pMEN.These results indicated that Proteobacteria were more sensitive to changes in marine variables, especially under high pCO 2 .They also indicated that no specific bacteria was significantly correlated with changes in either environmental variables or pCO 2 concentrations.

Discussion
Understanding the response of bacterioplankton communities to elevated pCO 2 is important for evaluating climate change effects on the ocean ecosystem.However, how the marine bacteria respond to elevated pCO 2 and decreased pH is not well understood and remains controversial (Liu et al., 2010;Joint et al., 2011).The EPOCA mesocosm experiments have provided a great opportunity to answer this critical question at the community-level (Riebesell et al., 2013b).The present study surveyed bacterioplankton communities by MiSeq sequencing of 16S rRNA gene amplicons and RMT-based ecological network analysis.Previous studies have demonstrated that phytoplankton respond significantly to elevated pCO 2 (Riebesell et al., 2013a, b;Schulz et al., 2013).Moreover, primary production and bacterial protein production are simultaneously stimulated by elevated pCO 2 (Engel et al., 2013;Motegi et al., 2013;Piontek et al., 2013), which generally results in net production decreases by the end of incubation experiments (Engel et al., 2013).Our results indicated that the phylogenetic diversity and structure of bacterioplankton communities were generally similar among different pCO 2 treatments, which is consistent with previous EPOCA bacterial community investigations (Roy et al., 2013;Sperling et al., 2013).Furthermore, the pMEN analysis suggested that the bacterioplankton community did not respond significantly to elevated pCO 2 at the community level, which means that our study agreed with the null hypothesis that changes in pCO 2 concentration have little influence on the affected microbial communities (Joint et al., 2011) and the bacterioplankton community had a certain resilience to the pCO 2 perturbations.The present study provides novel insights into how bacterioplankton communities respond to ongoing pCO 2 increases.
Generally, high diversity is thought to allow bacterial communities to resist disturbance, this is known as the "insurance effect" (Boles et al., 2004).Joint et al. (2011) suggested that high diversity is one of the mechanisms employed by bacterioplankton resist to elevated pCO 2 and decreased pH.Besides, micro-organisms do not live in isolation but form complex ecological interaction webs, which are also closely related to community stability (Faust and Raes, 2012).Changes global interactions (e.g.connectivity and modularity) among bacterioplankton communities along the pCO 2 gradient were unclear in the present study (Table 2 and Figure 2).These results suggested that the community interaction webs were relatively stable, despite the pCO 2 changes.In an ecological network, a negative correlation suggests a competitive relationship between species (Faust and Raes, 2012).The competition between species should yield less variation in total abundance or biomass, known as compensatory dynamics (Gonzalez and Loreau, 2009).A community with more competitors would be more stable under environmental fluctuations through maintaining relatively low covariations in community densities (Ives and Hughes, 2002).Klug et al. (2000) conducted a microcosm experiment and suggested strong compensatory dynamics when plankton were directly affected by pH perturbations.Furthermore, the competition might increase both the amplitude and asynchrony of population fluctuations generated by the asynchrony in species responses to the environment, which is the basis for functional compensation between species and potentially stabilizes aggregate community or ecosystem properties (Loreau, 2010).Ecosystem modelling indicates that asynchrony in species intrinsic rates of natural increase could stabilize the community (Fowler et al., 2012).Therefore, the dominant competitive relationships observed in our study (Figure 2) might sustain population diversity and enhance ecosystem resilience to perturbations.
Additionally, the high modularity observed in our study (Table 2) might result in a more complex ecosystem that would be more stable than dispersal ecosystems.Because of the tighter The nodes of Z i .2.5 and P i , 0.625 are indicated as the module hubs that were more closely connected within the module, while the nodes of Z i , 2.5 and P i .0.625 are the connectors that were more closely connected to nodes in other modules.Peripherals of Z i , 2.5 and P i , 0.625 are considered specialists in each module, while the network hubs of Z i .2.5 and P i .0.625 are super-generalists.
Bacterioplankton community resilience to OA interactions between nodes within the modules, the disturbances that influenced these nodes were unlikely to spread to other modules.For example, the OTU_1466 in the 175 matm pMEN only connected with nodes within that module, hence the abundance variation of OTU_1466 must affect more than two nodes to spread this variation to other modules (Figure 2).The extinction of a module hub might cause the module to fragment, having little or no cascading impact on other modules.In contrast, the extinction of connectors might cause the network to fragment but with minor impacts on the interactions within the isolated modules (Olesen et al., 2007).Therefore, compared with a non-modular structure, a modular structure should dampen the rapid spread of disturbance in a community (Olesen et al., 2007).Since each module in a pMEN indicates an ecological niche (Chaffron et al., 2010), our observation that negative relationships dominated between rather than within the modules suggested that the variations within a niche did not affect community stability.Therefore, competitive relationships and the modular structure of bacterioplankton communities provided resistance to variations in pCO 2 in the present mesocosm experiment.However, previous studies have shown that interactions among soil bacterial communities are influenced by elevated atmospheric CO 2 over a longer time-scale (i.e.annually; Zhou et al., 2010Zhou et al., , 2011;;Tu et al., 2015).Whether or not marine bacterioplankton also respond to long-term pCO 2 changes is worthy of investigation.
In an ecological network, the module hubs and connecters are considered keystone species, which have a crucial role in sustaining community stability (Paine, 1969).We found a relatively larger proportion of Gammaproteobacteria and Flavobacteria in pMEN module hubs and connectors in our study (Supplementary Figure S5 and Table S4).The relative abundance of Gammaproteobacteria and Flavobacteria did not change significantly among the pCO 2 treatments (Roy et al., 2013;Supplementary Figure S6) and they were not significantly correlated with geochemical variables under elevated pCO 2 .Therefore, these stable keystone species likely provided stability to the entire bacterial community.Monier et al. (2014) found the relative abundance of two groups of Gammaproteobacteria (Alteromonadales and Oceanospirillales) responded significantly to variations in pCO 2 in the north Arctic Ocean, but they did not respond significantly to decreased pH in terms of general phylogenetic structure.In our experiment, these potentially sensitive groups were not the major components and only a few OTUs belonging to these two Gammaproteobacteria groups were found in nine pMENs.Furthermore, we did not observe a clear trend in the interaction between the Gammaproteobacteria and their closest neighbours with elevated pCO 2 (Supplementary Figures S7 and S8).
Despite their general structure, bacterioplankton communities were insensitive to elevated pCO 2 in the EPOCA experiments, the low abundance species had the potential to respond to elevated pCO 2 .Roy et al. (2013) found that some rare species, like Methylotenera, Fluviicola and Colwellia, were significantly correlated with changes in pCO 2 .However, only a few of these bacterioplankton appeared within the nine pMENs in our analysis.For example, Colwellia appeared in the high pCO 2 treatment at the end of the incubation period, but was almost non-existent before that (Roy et al., 2013;Supplementary Figure S9).Because variations in keystone species have a greater impact on ecosystem stability and resilience than rare species (Paine, 1969), we used the major components of the bacterioplankton community to construct the networks.
Elevated pCO 2 had increased the POC and dissolved organic carbon (DOC) production by the end of the incubation period; therefore, bacterial protein activity and production were also stimulated (Motegi et al., 2013;Piontek et al., 2013).However, bacterial production, respiration, and growth efficiency were similar over the incubation period among the different pCO 2 treatments, which suggested that elevated pCO 2 had little impact on bacterial metabolism at the community level (Motegi et al., 2013).Previous meta-analyses have found no evidence that bacterial community growth efficiency is stimulated by elevated pCO 2 (Liu et al., 2010;Motegi et al., 2013).We did not find any significant correlation between network topology and bacterial production (except for the 675 matm treatment) and picophytoplankton abundance change.These results suggested that primary production stimulation had little influence on these keystone components.In fact, the variations in DOC and POC production during the incubation were greater those among the different pCO 2 levels (Engel et al., 2013).This is perhaps partially explained by the insignificant network response to pCO 2 among different mesocosms.However, another possible reason is different bacterial assemblage living strategies, i.e. free-living vs. particle-attached.Allgaier et al. (2008) reported that the free-living fraction of a bacterioplankton community was significantly correlated with elevated pCO 2 in a temperate fjord in Norway.However, the dynamics of the bacterioplankton community particle-attached fraction, compared with the free- living fraction, had a strong relationship with phytoplankton bloom, and exhibited a higher richness at the end of the incubation period in the EPOCA experiments (Sperling et al., 2013).The richness of the free-living fraction, however, decreased at the same time (Sperling et al., 2013).Thus, we would expect that the networks based upon entire bacterioplankton communities in the present study could not reflect the response of particle-attached bacteria to pCO 2 changes.Additionally, the lower bacterial abundance during this period combined with increased viral lysis in the higher pCO 2 mesocosm suggested top-down control on bacterial production (Brussaard et al., 2013), which might also not have reflected pMENs and suggested the importance of top-down control during this experiment.Therefore, the network analysis involved in both bacterioplankton and virioplankton, even zooplankton, would provide a more comprehensive explanation of the future effect of elevated pCO 2 in the ocean.
A few studies have found that decreased pH significantly influences specific bacterial groups in various oceans (Krause et al., 2012;Maas et al., 2013).However, their results were not consistent with ours; this might have been because of the difference between pH manipulations by HCl and CO 2 .Manipulation with HCl, which modifies the total alkalinity (TA) and pH with constant dissolved inorganic carbon (DIC), is not the most accurate approach to mimic future carbonate chemistry (Hurd et al., 2009;Krause et al., 2012).In contrast, CO 2 addition induces a pH decrease and DIC increase when TA is constant (Hurd et al., 2009).Even the carbonate chemistry parameters were similar between both manipulation methods when pCO 2 ,700 ppm (Hurd et al., 2009).However, when pH was modified to 7.5 by HCl, the [HCO − 3 ] value was 22% lower compared with the pH value obtained by pCO 2 manipulation (Hurd et al., 2009).The previous study has demonstrated that this difference in DIC affects some algal groups, e.g.those lacking a carbon-concentrating mechanism exhibited lower growth rate in pH 7.5 by HCl manipulation compared with pCO 2 manipulation (Hurd et al., 2009).Yet, the impact of this difference on the bacterioplankton is still unclear.Because mixotrophic micro-organisms are widespread carbon in the world's oceans (Hugler and Sievert, 2010), methodological differences may have affected the response of bacteria to changes in pH.Thus, pCO 2 manipulation is the most appropriate method available to mimic future changes in oceanic carbonate chemistry.

Conclusion
The hypothesis that "marine microbes possess the flexibility to accommodate pH change" is primarily based on the observation that microbial populations confront large variations in pH, both short-term and seasonally, in marine environments (Joint et al., 2011).This hypothesis is supported by recent studies at the microbial community level (Allgaier et al., 2008;Newbold et al., 2012;Roy et al., 2013;Sperling et al., 2013;Zhang et al., 2013).In combination with the present study, these data suggest that a larger population size and higher diversity within a microbial community contribute to the resistance of communities to high pCO 2 .Additionally, our network analysis provides further evidence that complex interactions (e.g.interspecies competition and modularity) among individual microbial species in a natural population help microbes resist and recover from pH and pCO 2 disturbances.This could explain the contradictory observations between OA effects on single species (e.g.pure culture) and natural marine microbe populations.
pCO 2 cosmosms contained M3 and M7, while elevated indicates the other 7 mecosoms.b Average phylum OTU percentage of total OTUs with standard deviation.c Average class OTU percentage of total phylum OTUs with standard deviation.868 Y. Wang et al.

Figure 1 .
Figure 1.Heatmap of sample compositions based on the proportion of each classes.Black represents the non-detected class; colour from dark to light indicates the proportion from low to high.The pCO 2 concentrations (matm) from A to I are 175, 180, 250, 340, 425, 600, 675, 860, and 1085.x-axis (D and a number) indicate the incubation day from 21 (DC-1) to day 30.

Figure 2 .
Figure 2. Bacterioplankton network interactions under 175 (a), 180 (b), 250 (c), 340 (d), 425 (e), 600 (f), 675 (g), 860 (h), and 1085 (i) matm pCO 2 .Five to seven modules were formed under different pCO 2 concentrations.Each node represents an out, which signified a species.Node colours indicate different phyla.Each line connects two OTUs.A blue line indicates a positive interaction between two individual nodes suggesting a mutualism or cooperation, while a red line indicates a negative interaction suggesting predation or competition.Top five highest-abundance OTUs are indicated by bigger dots and marked with OTU identification numbers.The cycle composed of several nodes is a module in the pMEN, which is more correlated in a module than between other modules.Only the module containing more than five OTUs are displayed.

Figure 3 .
Figure 3. pMEN submodules under nine different pCO 2 treatments.Each dot represents an OTU from nine different pCO 2 treatments (a -i indicate pCO 2 concentration from 175 to 1085 matm).The Z-P plot showing OTU distribution based on their module-based topological role according to within-module (Z) and among-module (P) connectivity.The nodes of Z i .2.5 and P i , 0.625 are indicated as the module hubs that were more closely connected within the module, while the nodes of Z i , 2.5 and P i .0.625 are the connectors that were more closely connected to nodes in other modules.Peripherals of Z i , 2.5 and P i , 0.625 are considered specialists in each module, while the network hubs of Z i .2.5 and P i .0.625 are super-generalists.
a Ambient indicates ambient

Table 2 .
Topological properties of the microbial communities' phylogenic molecular networks under nine pCO 2 concentrations and their rewired random networks.