Do invasive quagga mussels alter CO2 dynamics in the Laurentian Great Lakes?

The Laurentian Great Lakes have experienced unprecedented ecological and environmental changes, especially after the introduction of invasive quagga mussel (Dreissena rostriformis bugensis). While impacts on ecological functions have been widely recognized, the response of carbon dynamics to invasive species remains largely unknown. We report new CO2 data showing significant increases in pCO2 (up to 800 μatm in Lake Michigan) and CO2 emission fluxes in most of the Great Lakes compared to those prior to or during the early stage of the colonization of invasive quagga mussels. The increased CO2 supersaturation is most prominent in Lakes Huron and Michigan, followed by Lakes Ontario and Erie, but no evident change was observed in Lake Superior. This trend mirrors the infestation extent of invasive quagga mussels in the Great Lakes and is consistent with the decline in primary production and increase in water clarity observed pre- and post-Dreissena introduction, revealing a close linkage between invasive species and carbon dynamics. The Great Lakes have become a significant CO2 source to the atmosphere, emitting >7.7 ± 1.0 Tg-C annually, which is higher than the organic carbon burial rate in global inland-seas and attesting to the significant role of the Laurentian Great Lakes in regional/global CO2 budget and cycling.

Levels of carbon dioxide (CO 2 ) not only can serve as an indicator of autotrophic or heterotrophic nature of an aquatic environment 1 , but also can be an important parameter to elucidate calcification and potential pH changes in a water body [2][3][4] . Due to increasing human consumption of fossil fuels and other anthropogenic activities, atmospheric CO 2 level has increased from 280 μ atm in the 1800 s to levels above 400 μ atm in 2015 5 . As the most important greenhouse gas on earth, inventory and fluxes of CO 2 across different reservoirs and its role in the global carbon cycle and thus climate and environmental changes have been the focus of many recent research programs [6][7][8] .
Unlike the open ocean which contains abundant dissolved inorganic carbon (DIC, usually > 2000 μ mol/ kg) and typically serves as a sink of CO 2 7,9-11 , lakes and inland waters could serve as a source of CO 2 to the atmosphere [12][13][14] and may play an essential role to the local, regional and even global carbon cycles and climate. Recently, carbon cycling and magnitude of CO 2 fluxes across the air-water interface in lakes and inland waters have received increasing attention [15][16][17][18] , although studies on carbon dynamics in the Great Lakes remain scarce 17,19 .
The Laurentian Great Lakes are the largest freshwater system on Earth and receive vast amounts of organic and inorganic carbon from surrounding terrestrial ecosystems. Over the past decades, the Great Lakes have experienced significant ecological and environmental changes due to the introduction of invasive species, notoriously nonindigenous quagga mussels (Dreissena rostiformis bugensis), leading to a decrease in primary production and increasing water clarity in the Great Lakes [20][21][22][23][24][25] . In addition to impacts on foodweb structure and ecological functions, changes in biogeochemical cycles of nutrients have also been documented [26][27][28] . However, specific changes in carbon dynamics, as well as the impacts and biogeochemical consequences after the colonization of invasive quagga mussels in the Laurentian Great Lakes remain poorly understood. The direction and magnitude of air-lake CO 2 fluxes after dreissena introduction in the Great Lakes are largely unquantified, especially in Lakes Michigan, Huron and Ontario 17,[29][30][31] .
To examine the response of carbon dynamics to the introduction of invasive quagga mussels and linkages between invasive species and changes in biogeochemical cycling, open lake water samples were collected from all of the Laurentian Great Lakes, including Lake Superior, Lake Michigan, Lake Huron, Lake Erie and Lake Ontario, Scientific RepoRts | 6:39078 | DOI: 10.1038/srep39078 during August 2013 (Fig. 1). Furthermore, seasonal water samples were also collected from Lake Michigan between 2013 and 2015. In addition to water isotopic composition (δ 2 H and δ 18 O), total alkalinity (TA), dissolved inorganic carbon (DIC), and pH were measured prior to the evaluation of the partial pressure of CO 2 (pCO 2 ) and air-water CO 2 fluxes.

Results
Variations in pH, TA and DIC in the Great Lakes. The Laurentian Great Lakes can be characterized as a high-pH and high-carbonate ecosystem although Lake Superior had a relatively lower pH and carbonate abundance compared to the other Great Lakes ( Fig. 2 and Table 1). For example, pH values in Great Lakes waters were typically higher than 8 except for two sampling sites in Lake Superior (Table 1). Lakes Michigan, Erie and Ontario had higher pH values, averaging 8.22 ± 0.04, 8.19 ± 0.08 and 8.20 ± 0.04, respectively, followed by Lake Huron (8.07 ± 0.06) and Lake Superior (7.98 ± 0.05).
Concentrations of DIC and TA typically exceeded the threshold of 1000 μ mol/kg except for Lake Superior ( Table 1). Values of TA varied from 830 to 2200 μ mol/kg for all the Great Lakes, with an average of 843 ± 10 μ mol/ kg in Lake Superior, 2168 ± 22 μ mol/kg in Lake Michigan, 1600 ± 59 μ mol/kg in Lake Huron, 1872 ± 22 μ mol/ kg in Lake Erie and 1815 ± 19 μ mol/kg in Lake Ontario (Table 1). Similarly, different DIC concentrations were observed among the five Great Lakes, showing the highest value in Lake Michigan (average of 2065 ± 24 μ mol/kg) and lowest values in Lake Superior (average of 793 ± 26 μ mol/kg). Generally, these DIC and TA values are higher than those in most freshwater systems (e.g., refs 32, 33 and 34) and comparable to those of global ocean basins 35 (especially data of Lake Michigan, Table 1S). Our data here are also similar to the long-term data from previous studies in the Great Lakes 36 . High DIC and TA in Great Lakes waters are related to the spatial distribution of limestone and carbonate weathering in the Great Lakes basin. It is the high carbonate abundance that allows the thriving and rapid colonization of the mussel community in the Great Lakes 37 , which in turn has "re-engineered" the lake ecosystem during the past decades [22][23][24][25] .
The surface distribution ( Fig. 2) showed that TA, pH and DIC concentration all increased consistently from the upper Great Lakes (i.e., Lake Superior) to Lake Huron and then to the lower Great Lakes (i.e., Lakes Erie and Ontario), similar to the increasing δ 2 H and δ 18 O values in surface waters along the Great Lakes. This is consistent with the general water transport pathway from the upper to the lower Great Lakes, suggesting an accumulative effect on water chemistry in the lower Great Lakes from surrounding riverine inputs 38 . Among the five Great Lakes, Lake Michigan had notably high values for each of the measured parameters including pH, TA, DIC, and water isotopes ( Fig. 2 and Table 1), which is related to its high abundance of limestone and the fact that it is a semi-closed basin which allows the accumulative effect of evaporation in Lake Michigan.
Variations in pCO 2 along the Great Lakes. As shown in Table 1 and Fig. 3, our results consistently showed an evident CO 2 supersaturation (> 398 μ atm) for all the Great Lakes, with the highest pCO 2 in Lake Michigan (762 ± 88 μ atm) and Lake Huron (774 ± 92 μ atm), followed by Lakes Erie (725 ± 126 μ atm) and Ontario (647 ± 47 μ atm). Even in Lake Superior, with the lowest TA and DIC abundance (both < 1000 μ mol/kg) and lowest pCO 2 values (461 ± 77 μ atm), CO 2 supersaturation was still evident (Fig. 3). This indicated that surface waters of all the Laurentian Great Lakes can be a net CO 2 source to the atmosphere, consistent with those observed for other global lakes 16,18 .
Additionally, our seasonal pCO 2 data from Lake Michigan between 2013 and 2015 also showed a significant

Discussion
The invasive species, especially filter feeders including zebra and quagga mussels, have been shown to strongly affect the foodweb structure and biogeochemical cycles of nutrients after their introduction to the Great Lakes 21-25 . Thus, the equilibrium and interaction of CO 2 between lake and atmosphere may have been altered during the past decades, particularly after quagga mussels became predominant. As shown in Fig. 3, the pCO 2 values in the Great Lakes observed in 2013 are significantly higher than the median values of pCO 2 observed during the same sampling season between 1983-2006 17,30,39 (Table 1), especially for Lakes Michigan (p < 0.001), Huron (p < 0.001), Erie (p = 0.004) and Ontario (p = 0.035). Compared to the increase in atmospheric CO 2 levels from 1983 to 2013 (~60 μ atm), the increase in pCO 2 in the water column of the Great Lakes was considerably larger. Taking Lakes Michigan and Huron as examples, during 1983-2006, the average ∆ pCO 2 value, the difference in pCO 2 between atmosphere (averaging ~362 μ tm during that time periods) and lake waters, was about 118 and 22 μ atm for Lakes Michigan and Huron 30,39 , respectively. In comparison, ∆ pCO 2 values in these two lakes during summer 2013 reached a new high of ~400 μ atm, averaging 365 ± 89 and 377 ± 92 μ atm for Lake Michigan and Lake Huron, respectively (Table 1), showing an increase of 4 times in ∆ pCO 2 values in August between 1983-2006 and 2013 (Fig. 3). For Lakes Erie and Ontario, the increase in pCO 2 values between 1983-2006 and 2013 was also significant, although large errors existed in historical CO 2 data 30,39 for the shallower lakes with high primary production and anthropogenic influences including seasonal hypoxia and algal bloom, especially in Lake Erie 40,41 .
Moreover, our seasonal CO 2 data and its vertical distribution in Lake Michigan (Fig. 4) further demonstrated the increased CO 2 supersaturation in the water column. For example, pCO 2 in Lake Michigan during October 2015 ranged from ~600 μ atm in subsurface waters to ~770 μ atm in deeper waters (Fig. 4). Elevated pCO 2 in deeper waters resulted from the decomposition of organic matter and decreased biological uptake also implied a potential higher CO 2 emission flux from lake waters when the CO 2 -enriched deeper water is upwelled to surface during winter mixing seasons 42 . Furthermore, our seasonal data clearly show a consistent increase in pCO 2 from summer to fall and then to winter, reaching as high as 944 ± 83 μ atm during December 2015; an increase by ~30% compared to other seasons (Fig. 4). Thus, the pCO 2 data obtained during summer 2013 (Fig. 3) might represent the lower limit of CO 2 supersaturation in the Great Lakes, although other factors could have an influence, such as partial ice-cover during winter.
In contrast to Lake Michigan and other lower Great Lakes, the increase in pCO 2 in Lake Superior during summer between 1983-2006 17,30 and 2013 does not seem to be significant (p = 0.665, Fig. 3, Table 1). Nevertheless, increased CO 2 supersaturation in Lake Michigan and other lower Great Lakes coincides with the distribution and density of invasive quagga mussels among all the Great Lakes. For example, high quagga mussel densities (e.g., > 19,000-mussels/m 2 ) have been reported for Lake Michigan 22,25,27,43 , but a very low density for Lake Superior. Although zebra mussels arrived into the Great Lakes before quagga messuls, they have much weaker filter ability, and changes in foodweb structure and ecological function were more evident only after quagga mussels became predominant 25,[44][45][46] . In addition, quagga mussels did not start blanketing offshore regions of the Great Lakes, such as Lake Michigan until 2005 25,44 . Therefore, CO 2 data collected during 1983-2006 17,30,39 can be considered representative of the time period prior to or during the early stage of the colonization of invasive quagga mussels in the Great Lakes. Linking the spatiotemporal variations of quagga mussel population with the changes in the extent of CO 2 supersaturation during summer sampling time, we hypothesized that increasing CO 2 supersaturation in the Great Lakes during the past decade has been induced by the colonization of invasive quagga mussels and the subsequent biogeochemical response to that colonization.  Possible pathways and mechanisms to support our hypothesis include (1) decreasing primary production in the Great Lakes after the introduction of invasive quagga mussels 20-23,47,48 ; (2) increase in water clarity 24 which enables the light penetration into deeper waters and consequently enhances the photo-degradation of natural organic matter in the water column; and (3) metabolic processes of vast quantities of invasive quagga mussels blanketing lake floor 27,49 , although other processes, such as changes in sediment-water processes, coastal erosion and water chemistry after the colonization of quagga mussels, are also important. All these changes are inter-correlated and are the result of the colonization of invasive quagga mussels and all these processes would directly or indirectly favor the accumulation and release of CO 2 in the water column, resulting in increased CO 2 abundance in the Great Lakes, especially those heavily infested with quagga mussels (e.g., Lake Michigan 25 ). The close coupling between invasive species and carbon dynamics elucidated here in the Great Lakes clearly shows how small invasive quagga mussels could have caused basin-scale changes in carbon dynamics and biogeochemical cycling.
The increase in pCO 2 levels after the colonization of invasive quagga mussels in the Great Lakes evidently exceeded the increase in average atmospheric CO 2 level from 1983 to 2013 (~60 μ atm). Thus, enhanced CO 2 emission fluxes from lake waters to the atmosphere can be expected. Based on these summer pCO 2 data, the daily CO 2 emission fluxes ranged from 8.37 to 14.48 mmol-C/m 2 /d with an average of 10.38 ± 2.52 mmol-C/m 2 /d for Lake Michigan and ranged from 7.69 to 13.60 mmol-C/m 2 /d, averaging 10.71 ± 2.61 mmol-C/m 2 /d for Lake Huron (Fig. 5 and Table 1). These CO 2 emission fluxes were 3 to 4 times higher than those in Lakes Michigan and Huron at the same sampling time during 1983-2006, based on their differences in ∆ pCO 2 values. Lakes Erie and Ontario had CO 2 emission fluxes of 9.35 ± 3.59 mmol-C/m 2 /d and 7.13 ± 1.33 mmol-C/m 2 /d respectively, while an emission flux of 1.84 ± 2.17 mmol-C/m 2 /d was estimated for Lake Superior. Overall, the emission fluxes of CO 2 in all five Great Lakes follow the order of Huron ≥ Michigan > Erie > Ontario > Superior during summer 2013 ( Fig. 5 and Table 1), distinct from the order of Superior > Michigan > Huron > Ontario > Erie observed prior to or during the early stage of quagga mussel colonization 30,50 and further attests to the strong impact of invasive species on carbon cycles in the Laurentian Great Lakes.
Although CO 2 emission fluxes from the Great Lakes were comparable to or lower than those from small lakes (e.g., refs 51, 52 and 53), the large surface areas still made the Great Lakes a strong CO 2 source to the atmosphere. We estimated the lake-wide integrated CO 2 fluxes for each Great Lake and those estimated flux values are 1.5 × 10 8 mol-C/d for Lake Superior, 6.0 × 10 8 mol-C/d for Lake Michigan, 6.4 × 10 8 mol-C/d for Lake Huron, 2.4 × 10 8 mol-C/d for Lake Erie, and 1.4 × 10 8 mol-C/d for Lake Ontario (Fig. 5). Collectively, up to (7.7 ± 1.0) × 10 12 g-C/yr or 7.7 ± 1.0 g-C/yr in the form of CO 2 can be released to the atmosphere from the Laurentian Great Lakes, which is comparable to the annual DIC flux to the ocean from the Mississippi River and many other world rivers 54,55 . This annual CO 2 emission flux from the Laurentian Great Lakes is also higher than the annual organic carbon burial rate in all inland seas 16,56 , and makes up 1% of the annual CO 2 emission from global deforestation 57 . It should be noted that our estimated CO 2 emission flux could be the lower limit since it was calculated from summer sampling when pCO 2 levels were lower relative to winter seasons 42 (e.g., Fig. 4). Consequently, the Laurentian Great Lakes certainly serve as a significant CO 2 source to the atmosphere especially after the colonization of invasive species, and play an essential role in regional and global CO 2 budgets and cycling. This is a paradox in that the Great Lakes with similar carbonate abundance and pH as those in global oceans, should be expected to absorb CO 2 from the atmosphere and become acidified in the face of the . However, the Great Lakes have responded actively to the changing ecosystem and are making a positive feedback to climate change, with decreased whiting events (or CaCO 3 saturation/ precipitation) during summer 24 , increased CO 2 supersaturation and emission fluxes after the introduction of invasive quagga mussels and resultant decrease in primary production, enhanced degradation of natural organic matter due to increased water clarity, and direct CO 2 release from quagga mussel respiration 27,49 .

Materials and Methods
Sampling. During August 2013, surface water samples were collected from open lake stations on all of the Great Lakes: including Lake Superior, Lake Michigan, Lake Huron, Lake Erie and Lake Ontario (LS, LM, LH, LE and LO). Specific sample locations are shown in Fig. 1 and Table 1. Additionally, seasonal transect sampling, including a vertical profile at an open lake station (43°11.5362′ N; 87°39.9566′ W, 104 m depth), was also conducted in open Lake Michigan during 2013 to 2015. Detailed sampling locations can be found in ref. 28. Surface water samples were taken from a hull-mounted, all-Teflon and stainless steel, high-speed pumping system (2 m depth). Hydrographic parameters were recorded with a Hydrolab Datasonde 5A and coordinated with ship's positioning systems, such as temperature (°C) and chlorophyll-a fluorescence (V). Vertical profiles of hydrographic characteristics (SeaBird SBE 25Plus) and Niskin bottle-collected water samples were taken in open Lake Michigan at the offshore station. After collection, samples were saved for pH and total alkalinity (TA) and water isotopes,

Analysis.
Concentrations of DOC were measured on a Shimadzu TOC-L analyzer using the high temperature combustion method 59 . Total dissolved carbon (TDC) was measured on the same TOC analyzer without acidification and sparging. The concentrations of DIC were calculated by the difference between TDC and DOC concentrations 60 . Community consensus seawater reference from the University of Miami for DOC and the certified reference seawater (CRS) from the Scripps Institution of Oceanography for DIC 53 , as well as working standards were measured as samples to ensure data quality. After calibration with three standard solutions (pH = 4.00, 7.00 and 10.00), the pH electrode (Sartorius PB-11) was used to measure the pH. The precision and accuracy of pH was ± 0.01. Based on the Gran titration procedure 61  where k (cm/h) is the gas transfer velocity of CO 2 , and the K o (mol/m 3 /atm) is the solubility coefficient of CO 2 at in situ temperature and salinity. Values of k and K o were calculated based on the method of ref. 63 and ref. 64, respectively, and the monthly average wind speed in the study area was used for the calculation. The pCO 2-water and pCO 2-air are the partial pressures of CO 2 (in μ atm) in surface waters and air, respectively. The global annual averaged surface air pCO 2 for 2013 from NOAA (http://www.esrl.noaa.gov/gmd/) was used as the pCO 2-air value (398 μ atm). When the difference in partial pressure of CO 2 between water and air (∆ pCO 2 ) is higher than zero, the emission of CO 2 from surface waters to the atmosphere occurs. In contrast, a negative value of ∆ pCO 2 indicates uptake of atmospheric CO 2 in surface waters.
Uncertainties of pCO 2 and statistical analysis. Both DIC-pH and TA-pH pairs were used to derive the pCO 2 values in lake waters. Uncertainties of pCO 2 derived from pH ( ± 0.01) were estimated to be ± 16 μ atm. For the discussion on CO 2 fluxes, pCO 2 data derived from pH-DIC were used, while data derived from pH-TA were used to evaluate possible contribution from non-carbonate alkalinity 65 . As shown in Table 1, pCO 2 values calculated from pH-TA pair were, on average, 48.8 ± 14.4 μ atm higher than those derived from pH-DIC pair.
In other words, overestimation of pCO 2 from non-carbonate alkalinity could be up to 6.5 ± 1.4% in Great Lake waters even under the generally low DOC concentrations (Fig. 2). A one-way ANOVA analysis was performed to determine the significance of differences in pCO 2 between different sample groups.