Effects of acidification on the proteome during early development of Babylonia areolata

Increases in atmospheric CO 2 partial pressure have lowered seawater pH in marine ecosystems, a process called ocean acidification (OA). The effects of OA during the critical stages of larval development may have disastrous consequences for some marine species, including Babylonia areolata (Link 1807), a commercially important sea snail in China and South East Asia. To investigate how OA affects the proteome of Babylonia areolata, here we used label‐free proteomics to study protein changes in response to acidified (pH 7.6) or ambient seawater (pH 8.1) during three larvae developmental stages of B. areolata, namely, the veliger larvae before attachment (E1), veliger larvae after attachment (E2), and carnivorous juvenile snail (E3). In total, we identified 720 proteins. This result suggested that acidification seriously affects late veliger stage after attachment (E2). Further examination of the roles of differentially expressed proteins, which include glutaredoxin, heat‐shock protein 70, thioredoxin, catalase, cytochrome‐c‐oxidase, peroxiredoxin 6, troponin T, CaM kinase II alpha, proteasome subunit N3 and cathepsin L, will be important for understanding the molecular mechanisms underlying pH reduction.

Increases in atmospheric CO 2 partial pressure have lowered seawater pH in marine ecosystems, a process called ocean acidification (OA). The effects of OA during the critical stages of larval development may have disastrous consequences for some marine species, including Babylonia areolata (Link 1807), a commercially important sea snail in China and South East Asia. To investigate how OA affects the proteome of Babylonia areolata, here we used label-free proteomics to study protein changes in response to acidified (pH 7.6) or ambient seawater (pH 8.1) during three larvae developmental stages of B. areolata, namely, the veliger larvae before attachment (E1), veliger larvae after attachment (E2), and carnivorous juvenile snail (E3). In total, we identified 720 proteins. This result suggested that acidification seriously affects late veliger stage after attachment (E2). Further examination of the roles of differentially expressed proteins, which include glutaredoxin, heat-shock protein 70, thioredoxin, catalase, cytochrome-c-oxidase, peroxiredoxin 6, troponin T, CaM kinase II alpha, proteasome subunit N3 and cathepsin L, will be important for understanding the molecular mechanisms underlying pH reduction.
The gradually increasing atmospheric CO 2 partial pressure (pCO 2 ) has lowered the seawater pH in marine ecosystems. This change in seawater pH is regarded as ocean acidification (OA). The pH of the ocean's surface waters is predicted to decrease by 0.3-0.4 units by the year 2100 [1]. An estimated reduction of up to pH 7.3 in sea water is expected within the next 300 years [1][2][3].
Marine shellfish is a key community that stabilizes the marine ecosystem, is considered as the driving force of the ocean material flow and energy flow, and plays an important role in the marine carbon cycle [4]. The monetary value of OA impact on overall shellfish production is estimated to be about 100 billion USD [5]. The effects of a changing living environment on larvae are important, especially because these early development stages often represent a bottleneck for species struggling to survive in the changing marine environment [6,7]. The effects of OA during the critical stages of larval development may have disastrous consequences for some marine species [8], for instance, inhibition of the fertilization, embryonic and larval development, and survival in the acidified seawater of Sakura clam (Macoma balthica) [9]. Gazeau et al. [10] found that OA reduces the hatching and adhesion rates of mussels (M. edulis). OA can delay early embryo and larval development and show teratogenic and lethal effects in many marine invertebrates [11], including California mussels (Mytilus californianus) and Ruditapes decussatus [12,13]. The reductions of 0.7 units often lead to an abnormal development and high mortality in some shellfish, for instance, C. gigas [11], Mytilus galloprovincialis [14,15], and Pinctada martensii [16].
Differences in growth and larval development time may have important consequences for population dynamics and demography [17]; larval developmental plasticity is a crucial source of variation in adult phenotypes and can directly influence the evolution of populations and species [18,19]. Larval development includes a pelagic larval phase and benthic phase in some gastropod; larvae undergo metamorphosis and transited to benthonic larval phase and juveniles [20]. The attachment and metamorphosis of gastropod larvae affect the quantity change, population distribution, and species evolution [21]. Therefore, investigation on changes of gastropod larval development in pH may provide a realistic indication of how shellfish respond to OA.
Overall, OA can negatively affect the main physiological processes of shellfish, OA can influence calcification, destroy the acid-base balance of organisms, change metabolism, and reduce the immune function. However, the current related research methods are mainly based on traditional physiology and ecology [22], but lack studies which focus on the molecular response mechanism. Using molecular '-omics' techniques to reveal the molecular signals underlying the response of development characteristic to changes in pH is now possible [23]. Proteomics can be a comprehensive and accurate method to analyze the response mechanism of marine organisms to OA, the fitness of an organism can be assessed by proteome-molecular traits, which are directly related to their physiological performance [24]. Therefore, identification and characterization of proteins involved in various metabolic processes are required to predict the capacity of the organism to adapt to OA.
In the past decade, the environmental proteomics studies lead to the identification of molecular pathways and differentially expressed proteins of shellfish in response to OA [25,26]. Proteomic studies have shown that elevated pCO 2 (resulting in pH 7.3-7.87) leads to the differentially abundant proteins associated with oxidative stress, protein synthesis, energy metabolism, and the cytoskeleton in the oysters Crassostrea virginica, C. hongkongensis, C. gigas, and Saccostrea glomerata [24][25][26][27][28][29][30][31][32][33], and affects energy and primary metabolisms, stress responses, calcium homeostasis, the nucleotide metabolism, and cytoskeleton structure in pacific oyster Crassostrea gigas [34]. Identifying these proteins as well as determining the effect that OA exerts on oysters is fundamental to understand the mechanisms involved [25]. Dineshram et al. [29] examined the effects of OA on the proteome of the oyster species Crassostrea hongkongensis; over 7% of the proteome was altered in response to OA at pH 7.6. Proteomic analysis of wild-type C. gigas and B2 S. glomerata (a selective breeding line) showed that most of identified proteins were upregulated in the wild-type population. However, these proteins were downregulated in the B2 S. glomerata by CO 2 stress (pH 7.5-7.84) [26]. OA could reduce the expression level of some proteins to conserve energy, but OA could induce oxidative stress proteins [35]. Some of these OA-responsive molecules and pathways have been identified in marine invertebrate organisms, including oysters [24,25,36].
In proteomic of pCO 2 on larvae development, proteomics was used to compare the global protein expression pattern of oyster larvae exposed to ambient and to high CO 2 [27]. Dineshram et al. [28] proteomic studies revealed a set of proteins, referred to as 'protein expression signature', that changed their expression when the larvae were exposed to OA. Comparative proteome analysis data confirmed that there was a significant suppression of molecular pathways involved in protein degradation, energy production, and tissue growth, which correlate with decreased larval metamorphic behavior after the pH stress at 7.4, and more than 11% proteins of larval were downregulated with decreased pH [30]. Similarly, a study about the proteome in an earlier larval stage of C. gigas showed that a third of expressed proteins were downregulated in response to decrease pH [27]. Larval forms of tubeworms, barnacles, and oysters also showed increased downregulation of proteins associated with calcification, energy production, and cytoskeleton in response to decreased pH (pH 7.6) [37,38]. However, the specific function of these proteins in larva as a tolerance molecular mechanism to decreased pH is unknown.
The ivory shell, B. areolata, is the most promising economic marine gastropod in the present century due to their wide distribution along most parts of Asian coastlines and the southeast coast of mainland China, their delicate flavor, and high market reception. The annual output of B. areolata is more than 1000 tons, and the worth is more than 100 million RMB Yuan [39,40].
In this study, three different stages of larvae of B. areolata were cultured at high (OA) and ambient CO 2 (control) conditions, larvae of similar physiological age and size from the OA and the control groups were analyzed by label-free quantitative proteomics, and differentially expressed proteins were identified by mass spectrometry (MS). We investigated OA on the larvae of B. areolata as a model to address the following questions. How does the larval proteome respond to OA? Which proteins are involved in the response? Our objectives were as follows: This labelfree quantitative proteomics approach will enable us to partially elucidate as to how calcification, metabolism, oxidative stress, and stress tolerance proteins were modified in response to OA. Furthermore, to our knowledge, there is no other study to date using proteomics to differentiate OA-exposed B. areolata larval proteome versus healthy controls. Proteomics can identify crucial proteins involved in larval developmental phenotypes in response to future changes in pH and provides insights into the potential mechanisms involved in stress response to and tolerance of high CO 2 levels.

Experimental apparatus and organisms
The construction of seawater CO 2 system was advocated as previously described [41], which was the study result of our laboratory. The pCO 2 manipulation system is shown in Fig. S1. The methods are briefly described as follows: The principle of experimental system was to mix CO 2 -free, dry air, and pure CO 2 (99.99% purity) together at different ratios using mass flow rates to produce CO 2 -enriched air with different pCO 2 [42]. Atmospheric air was provided from an oil-free, medical air compressor. The water and particles were removed through two filters (GFR600-25, AIRTAC), then two polypropylene (PP) columns filled with soda lime to absorb the CO 2 were connected, and then a similar PP column filled with anhydrous CaCl 2 to further remove water was connected again. Finally, a 5-cm disk-type air filter was connected before connecting the regulating utilities. The treated air was then delivered into a pressure regulation valve and then into a needle valve, which can maintain a stable air flow.
Mass flow of CO 2 and air was measured by two mass flow sensors and was set to the desired level of CO 2 results in a CO 2 concentration of 800 ppm by adjusting the needle valves. The pCO 2 manipulation and measurement system were performed as previously described in the literature (see our previous study [41]). The method is described as follows (see Fig. S1): both CO 2 and air were homogenized at the bottom of a plastic container to generate CO 2 -enriched air with different CO 2 concentrations. A small proportion of these gas mixtures were directly monitored with a CO 2 detector (Li-7000, LI-COR) through a bypass.
Further adjustment of the needle valve was carried out before the reading system of the detector achieved the preset standard. The mass flow sensors and CO 2 detector were connected to a computer, and measurements were recorded with the associated software. The majority of CO 2 -enriched air was mixed with filtered seawater (0.22 lm filtered). Two small holes were drilled in the lid of the plastic container for input and outlet of CO 2 -enriched air. The container was kept semiclosed in order to ensure that the inner space of the container was filled with CO 2 -enriched air and the used seawater which was separated from the ambient air. Carbonate system parameters were calculated using CO 2 SYS software (version 2.1, Oak Ridge, TN, USA). Temperature and pH were measured with a pH meter. A total of 2 CO 2 concentrations were used: 400 ppm, which represents the current CO 2 concentration, and 2000 ppm [43], which represents the predicted CO 2 concentrations for the year 2300.

Samples acquisition
Adult B. areolata broodstock were purchased from the Haicang farm, Fujian province, and individuals with vitality were selected. Snails were maintained in an ecological culture pond (27.3 AE 1.1°C, salinity 30.1 AE 0.9, pH 8.0 AE 0.65, oxygenated water). The snails were fed daily with oysters and chopped fresh fish. All experiments involving animals reported in this study were approved by the Animal Care and Use Ethics Committee of the Xiamen University.
After B. areolata broodstock spawning, the oocysts were washed with filtered seawater and then transferred to 40-mL sample bottles containing filtered seawater with different pH values (pH 7.6 and 8.1). One oocyst was placed in each sample bottle and changed sea water one time each day until the larvae hatched from the oocysts. In order to ensure the stable pH value of the water, the sample bottle was put into the preservation box and CO 2 -enriched air was input into the box. The incubation time was recorded, and the samples were collected after the larvae were hatched from the oocyst. Under continuous exposure, embryos were cultured at different CO 2 concentrations from fertilization until the end of the experiment. For B. areolata, fertilization, hatching, and developmental conditions were measured.
In this study, three developmental stages of B. areolata, the veliger stage of B. areolata before attachment, the late veliger stage of B. areolata after attachment, and the juvenile stage (see our previous study [40]) were sampled. When larva is the veliger larval prior to settlement stage, many changes ensue, including the loss of velum and gill development. Then, with the development of foot, the larvae crawled and began their feeding. When morphological changes occurred, with the transition stage from planktonic larvae stage to benthic stage, that was from the late veliger stage to the juvenile stage, the feeding habits of B. areolata also gradually changed from phytophagy to sarcophagy. The larvae of every developmental stages were divided into five replicates. The control group, C, was at 400 ppm, and pH was at 8.15 AE 0.01, whereas the experimental group, E, was at 2000 ppm and pH at 7.61 AE 0.02. The water temperature was kept at 27°C. For the veliger stage before attachment, the control and acidified group were named C1 and E1, respectively; for the late veliger stage of metamorphosis at the bottom (velum atrophy), the control and acidified group were named C2 and E2, respectively; for the carnivorous juvenile snail, the control and acidified group were named C3 and E3, respectively; three larvae stages are shown in Fig. S2A-C. Early veliger larvae were fed daily with algae. When the late veliger stage sinks to the substratum and begins a benthic life, the feeding habits of larvae also gradually changed from phytophagy to sarcophagy; then, the snails were fed daily with oysters. Under these culture conditions, >85% of embryos developed into early veliger stage in 60 h, after which, the seawater was replaced with new seawater, larvae were cultured for 2 days, and then the veliger larva enters the late veliger larva stage. After 11-12 d, the late veliger stage enters the juvenile snail stage. Larval samples with similar physiological age and shell size between the control and the OA treatment were used for comparative proteome analysis. For every development stage, five replicates were set at each pCO 2 level. Three biological duplicates (triplicates) were selected from each of the five replicates for proteome analysis. Samples collected three development stages fulfilled with those similar larval size and age requirement for proteome analysis.
In every development stage, almost all individuals were in the same development period by microscopic examination. Approximately 4000-5000 of larvae per stage, fresh and living larvae, were filtered using an 80-lm filter, and then, the larvae were collected into 1.5-mL tubes for further analysis. The samples were centrifuged at 13 000 g for 3 min at 4°C. The larval material (15 lL of pellets) was stored at À80°C. For every development stage, the larvae were divided into three independent samples, and three replicates were performed for each development stage to ensure reproducibility.

Total protein extraction with TCA/acetone
A total of 100 mg sample was suspended immediately in 1 mL 10% TCA/acetone (w/v) solution, mixed, and then stably placed over 2 h at À20°C. The sample was centrifuged at 17 000 g for 15 min at 4°C; the supernatant was removed and then added with 1 mL 10% TCA/acetone (w/v) solution with 20 mM DTT. The sample was solubilized with a sonicator (BRANSON, S-450D), the sample was placed on ice, the program of sonicator was run by using 5-10 min blasts of 15% amplitude lasting 5 s each, with 4 s breaks, and then stably placed at 4°C for 30 min. Samples were fully homogenized with a tissue homogenizer and stably placed at 4°C for 30 min. The sample was centrifuged at 17 000 g for 15 min at 4°C, and the supernatant was removed. Then, the protein pellet was washed three times using 80% acetone solution with 20 mM DTT and added 1 mL each time. The sample was placed at 4°C for 30 min. Then, the sample was centrifuged at 17 000 g for 15 min at 4°C; the upper layer was removed, and the lower layer pellet was sampled. The precipitation block was carefully agitated with sterile toothpick every time. The protein pellet was washed with 1 mL 20 mM DTT 100% acetone and then placed at 4°C for 30 min. The sample was centrifuged at 17 000 g for 15 min at 4°C; then, the supernatant was removed, and the protein pellet was sampled. Then, the protein pellet concentration was redissolved in solution (6 M urea, 2 M thiourea, 0.1% SDS, 1 mM DTT, 10 mM Hepes, pH 8.0). The concentration of protein was measured using a Protein 2-D Quant kit (GE Healthcare, Piscataway, NJ, USA).

Digestion of proteins
In brief, 50 lg of proteins was alkylated with 55 mM iodoacetamide for 30 min at room temperature in the dark. Then, proteins were incubated with trypsin (1 g/50 g protein) (Promega, Madison, WI, USA) for 12 h at room temperature. The digestion of proteins was stopped by adding 5% formic acid (the final concentration was about 0.1%).

Desalting with C18 tips
Desalination was done using prewet tips with 100 lL of 50% acetonitrile (ACN) solution; the tips were discarded after. This step was repeated one time. The tips were then balanced with 100 lL of 0.1% trifluoroacetic acid (TFA) and repeated once. The digestion sample of 100 lL was sucked using the tips and then blown out; this step was repeated more than 10 times. Finally, the samples were blown out, and then, 5% acetonitrile containing 0.1% TFA was sucked with the tips, repeated once and followed elution in 100 lL of 50% acetonitrile containing 0.1% formic acid. The collected eluate was washed two times every time; 100 lL 95% acetonitrile containing 0.1% formic acid was sucked. After concentration and drying of elution liquid, the sample was then fully dissolved with 30 lL liquid (2% methanol and 0.1% formic acid) by vortex; polypeptide concentration was measured with Nanodrop microspectrophotometer, stored at À70°C, or centrifuged at 14 000 g for 15 min, and then directly analyzed.

Strong cation exchange liquid chromatography
The first dimension was separated using SCX column; the flow rate was 1000 lLÁmin À1 , divided into 10 compositions according to the chromatographic peak, then desalted by C18 Cartridge and lyophilization, and then frozen at À80°C. Gradient of liquid chromatography is set as follows: 0-25 min, 100% buffer A (KH 2 PO 4

Nanoacquity ultraperformance liquid chromatography
The two-dimensional sample was separated using Easy nLC liquid-phase system. The sample was dissolved with 50 lL solution (0.1% formic acid, 5% acetonitrile) and injected with autosampler into the Zorbax 300SB-C18 peptide traps (Agilent Technologies, Wilmington, DE, USA). Then, the sample was separated using chromatographic column; the column flow rate is set to 200 nLÁmin À1 . Solution A (0.1% formic acid aqueous solution) and solution B (0.1% formic acid acetonitrile water solution, in which acetonitrile was 100%) were used as the mobile phases with the gradient elution procedures. The liquid-phase gradient was as follows: 0-100 min, the range of B linear gradient was 5%-30%; 100-130 min, the range of B linear gradient was 30%-40%; 130-135 min, the range of B linear gradient was 40%-90%; and 135-140 min, the range of B was maintained at 90%. Chromatographic column was 75 lm 9 150 mm (RP-C 18 ), which was equilibrated with 95% solution A.

Q-Exactive mass spectrometric protein identification
Raw MS/MS data files were processed to peak lists in PRO-TEOME DISCOVERER software (Thermo Fisher Scientific, software version 1.4, Waltham, MA, USA) and then submitted to Mascot-generated files searches (Matrix Science, London, UK). The data-dependent tandem mass spectrometry labelfree analysis was carried out on a Q-Exactive mass spectrometer (Thermo Fisher Scientific, MA, USA); the flow rate was 300 nLÁmin À1 . Full scan spectra were acquired with a maximum ionization time of 200 ms, the spectra were recorded at 70 000 resolution, and the AGC target value was set as 10 6 ions; the number of scan ranges was 1, and then, 10 most abundant peaks for MS/MS were selected using a 18 000 resolution scan (the AGC target value was set as 1 9 10 4 ions; the maximum ionization time was 200 ms) with a 1.5 m/z ion selection window, and a normalized collision energy of 30% was used for collision-induced dissociation. A 40 s dynamic exclusion window was used to avoid repeated selection of peptides for MS/MS. MS2 activation type: HCD; MicroScan: 1; Underfill ratio: 0.1%.

Analysis of the Q-Exactive mass spectrometric data
Thermo Fisher Scientific Orbitrap RAW files were converted into complete peak lists, RAW spectra were converted into Mascot-generated files (mgf) using PROTEOME DISCOVERER software (Thermo Fisher Scientific, software version 1.4) and analyzed by using SIEVE software (Thermo Fisher Scientific) to quantify all detected peaks. All data were acquired based on 99% confidence for protein identification by the false discovery rate (FDR) ≤ 1%. Statistical analysis was carried out using a one-way ANOVA. P-values ≤ 0.05 acquired by Tukey's test were considered significant. The q-value was used to estimate false-positive results. Peaks exhibited expression changes (>1.5-fold). The search parameters were as follows: the enzyme specificity was trypsin; database: gastropods; the allowance was two missed cleavage site; the variable modification was Met oxidation; the fixed modification was carbamidomethyl (cysteine); protein mass was unrestricted; monoisotopic mass values were obtained; the fragment mass tolerance was AE 0.2 Da; and the peptide mass tolerance was AE 30 ppm. Protein FDR ≤ 0.05.

Bioinformatics analysis
The Gene Ontology (GO) database provides a controlled vocabulary for describing genes and gene products, as well as provides full access to this information in several formats. UniProt Knowledgebase (UniProtKB) was used for protein function analysis and acts as the central hub of protein knowledge by providing a unified view of protein sequence and function. Currently, UniProtKB database included more than 78 037 000 terms. The web site is http://www.uniprot.org/uniprot/.

Principal component analysis
Principal component analysis (PCA) is a statistical procedure that converts a set of observations of possibly correlated variables. The differentially expressed proteins characteristic parameters in six groups were analyzed by PCA (Fig. 1); principal component 1 (PC1) and principal component 2 (PC2) were obtained. The load value of six groups characteristics to two principal components and the variance contribution proportion of two principal components are shown in Fig. 1A,B. In Fig. 1A, the contribution proportion of PC1 was 43.3%, the contribution proportion of PC2 was 15.0%, and the PC1 and PC2 cumulative contribution proportion was 58.3%, biological duplicate samples were clustered together, the results showed that there was little difference between intragroups, and there was good repeatability between biological duplication samples. The results were extremely variable between C2 and E2 (Fig. 1B), suggested that C2 larva stage was sensitivity to OA, and acidification seriously affects late veliger larvae after attachment.

Identification of differentially expressed proteins
In this study, a total of 720 proteins (FDR < 1%, at least two unique peptides) were identified, and 245 differential abundance proteins (>1.5-fold) were acquired. As shown in Table 1, the protein abundance of E group is regarded as upregulation or downregulation. The results showed that the total number of differentially expressed proteins was 148 between C2 and E2 groups, whereas the number of proteins was 57 between C1 and E1 and was 26 between C3 and E3. Differentially expressed proteins were determined using normalized % protein amount, with protein identifications summarized in Tables S1-S3. For C1 and E1, eight differentially expressed proteins that increased over five times in E1, such as histone 2A (12.5-fold), glutaredoxin (7.14-fold), tubulin beta chain (6.25-fold), cathepsin D (5.56-fold), huntingtin interacting protein K (5.56-fold), catalase (5.56fold), and cytochrome-c-oxidase subunit I (5.26-fold). Seven differentially expressed proteins that decreased over three times in E1 included tumor rejection antigen-like protein (17.37-fold), 26S protease regulatory subunit 6B (8.36-fold), pedal peptide precursor protein  (Tables S2). For C3 and E3, nine differentially expressed proteins that increased over 2.20 times in E3 included tumor rejection antigen-like protein ( (Table S3). For C1 and E1, identification of differentially expressed proteins involved in the biological process is summarized in Fig. 2. The upregulated proteins in E1 were involved in response to oxidative stress, transduction and head involution. The downregulated proteins in E1 were involved in protein polyubiquitination,  detection of calcium ion, carbohydrate metabolic process, and ATP binding. For C2 and E2, identification of differentially expressed proteins involved in the biological process is summarized in Table S4. The upregulated proteins in E1 were involved in response to oxidative stress, head involution, protein glycosylation, protein dephosphorylation, lipid metabolism, and ATP synthesis; the downregulated proteins in E2 were involved in stress response and immunity.
For C3 and E3, identification of differentially expressed proteins involved in the biological process is summarized in Fig. 3. The upregulated proteins in E3 were involved in protein biosynthesis, NADH oxidation, and transduction; the downregulated proteins in E3 were involved in signal translation and tricarboxylic acid cycle.

Series protein expression data analysis by the Short Time-series Expression Miner
Short time-series expression miner is a tool for the analysis of short time-series gene (or protein) expression data [44], which can also analyze protein expression. The data from the B. areolata larva sampled at six groups (C1, C2, C3, E1, E2, and E3) were filtered to contain the 245 differentially expressed proteins that exhibited 1.5-fold upregulation or downregulation for at least one sample; the number of types was set 20. Comparing with the C1, all samples relative to C1 were calculated the log2 value.
Along each window are statistics on the number of proteins which were assigned to the profile and the enrichment p-value. The number at top of a profile box shows the profile ID number. The colored profiles had a statistically significant number of proteins assigned ('genes' in the graph refer to 'proteins' data). Colored profiles of the same color represent profiles grouped into a single cluster (Fig. 4). The bottom of each colored window shows profiles with P-value enrichment for proteins annotated (P-value < 0.05). Differentially expressed proteins in profiles 3, 7, 12, 13 and profiles 17 were significant (Fig. 5). Differentially expressed proteins in profiles 3, 7, 12, 13 and profiles 17 were summarized in Tables S5-S9. For profile 3, proteins with differential expression included ATP synthase beta subunit, putative mitochondrial ATP synthase F chain, mitochondrial malate dehydrogenase precursor, gelsolin, endo-1, 3-beta-D-glucanase, p38 MAPK, actin, putative tubulin beta chain, tektin A1, axonemal dynein light-chain p33, collagen pro alpha-chain, munc18-1-interacting protein 1, poly-(ADP-ribose) polymerase I, and non-neuronal intermediate filament protein A.

Hierarchical cluster analysis
Hierarchical cluster analysis was used to examine the patterns of similarities and differences among genotypes. Protein abundant hierarchical clustering was performed to identify main classes of protein abundance. Distances were determined for the six groups (C1, C2, C3, E1, E2, and E3) with a TMEV heatmap, and a dendrogram was constructed using HEMI 1.0 statistical software. Distances were calculated based on protein value intensities in all technical replicates. In this study, 720 proteins were detected; we selected a mean value-based approach and then proceeded with analysis. Results showed that C2 and C3 clustered in one clade. The result showed that formed clusters indicated clear separation of E2 from other groups (Fig. 6).

Discussion
Label-free quantitative proteomics provided rigorous and powerful tools for analyzing protein changes. In the present study, 720 proteins were identified. The low rate of solid attribution may have been caused by the limited cDNA information available for B. areolata or by improper annotations of B. areolata genes. B. areolata is a nonmodel species; the tryptic digestion resulted in acceptable MS/MS spectra; however, these mass signals did not match any peptide available [45] in the B. areolata or other gastropods identifiable protein databases. In this study, the number of differentially expressed proteins between C2 and E2 was largest; this finding suggests that acidification seriously affects late veliger larvae after attachment.
Effects of OA on proteins associated with respiration and oxidative stress In this study, the expression level of differential proteins associated with oxidative stress increased in E1 and E2; NADH is an important material for cell respiratory metabolism, and NADH oxidation function was upregulated in E3. The result showed that the number of upregulation proteins which are located on the mitochondria is larger than that of downregulation proteins in E1, E2, and E3.
The previous study showed that the metabolic process (including food intake, oxygen consumption, and ammonia excretion) of three bivalves: the pearl oyster Pinctada fucata, the noble scallop Chlamys nobilis, and the green-lipped mussel Perna viridis was likely affected by the reduction in seawater pH [46]. Fern andez-Reiriz et al. [47] found that the acidification conditions tested significantly reduced the respiration rates and increased the ammonia excretion rate ofRuditapes decussatus seeds. Zhang et al. [48] found that respiration rate of Chlamys farreri also decreased significantly at pH 7.3. Ding et al. [49] found that the respiratory metabolism pathway and the concerning enzymes of marine organisms might be changed as ocean acidification development. Those studies suggest that marine organisms have a series of adaptive measures to rescue the negative effects of low pH on respiration; these measures include increasing breathing frequency, expanding the depth of breathing, oxygen consumption increased first and then decreased, and a series of complex physiological regulation mechanism.
In the present study, cytochrome-c-oxidase subunit I increased 5.26-fold in E1; cytochrome-c-oxidase is the terminal enzyme in the electron transport chain reaction. Thompson et al. [26] showed that proteins involved in the mitochondrial electron transport chain and oxidative stress were the most affected by OA. In the identification of proteins between E1 and C1, the amount of glutaredoxin increased 7.14-fold among the upregulated proteins in E1, and catalase increased 5.56-folds in E1. Between E2 and C2, thioredoxin increased 85.70-fold in E2. Glutaredoxin and thioredoxin are considered to be a typical antioxidant, and thioredoxin can also remove the hydrogen peroxide; the catalase catalyzes the decomposition of hydrogen peroxide into oxygen and water.
Matoo et al. [50] studies have already demonstrated the negative impacts of elevated temperature and CO 2 levels on the metabolism and cell oxidation in these bivalve species. The combined effects of decreased pH and increased temperature values on antioxidant enzyme activities and lipid peroxidation were evaluated in the clam Chamelea gallina and the mussel Mytilus galloprovincialis [51]. Superoxide dismutase, catalase, and glutathione S-transferase activities as well as lipid peroxidation were measured. The results demonstrated that the experimental conditions significantly influenced the biochemical parameters of the bivalves [51].
In this study, peroxiredoxin 6 increased 2.04-fold in E1 (Table S1), peroxiredoxins (Prxs) work as a redox sensor to resist peroxide attacks. Previous studies showed that upregulation of some Prxs can protect cells from oxidative damage [52]. A similar study also showed the significant variation of Prxs in mantle of eastern oyster C. virginica and hepatopancreas of pacific oyster C. gigas under elevated pCO 2 condition using proteomic analysis [34,35]. Therefore, the upregulation of Prx 6 perhaps suggested the stress responses to elevated pCO 2 exposure, especially oxidative stress.
Effects of OA on proteins associated with response to immunity Two differential expression proteins were identified as heat-shock protein 70, which decreased 23.13-fold and 5.50-fold in E2, respectively (Table S2). Dineshram et al. [30] showed that decreased pH caused the reduction of several heat-shock proteins (Hsps), and the downregulation of Hsps has also been documented in barnacles under decreased pH [37] and sea urchins [53].
In GO of the biological process involved the identification of proteins. In comparison with C2, the protein involved in the stress response and immune function in E2 decreased. Organisms in the early developmental stages, namely, larval and juvenile, have low immunity that may be affected by the degree of acidification [54]. Gazeau et al. [10] and Matozz et al. [51,55] showed that acidification can significantly reduce the immunity in Mytilusedulis, Chamelea gallina, and M. galloprovincialis and decrease the rate of metabolism. Michaelidis et al. [56] showed that water pH below 7.5 was harmful on the immune function of mollusk; this result leads to a permanent reduction of hemolymph pH value, growth, and metabolism; however, the level of hemolymph bicarbonate from shell dissolution increased. Berge et al. [57] showed that acidification made a negative effect on the immunity and growth of shellfish; they also suggested that the decline of immune function was related to the negative metabolism. Prolonged exposure to pH reduction of 0.5-0.7 pH units possibly elicits strong stress responses in many organisms [15].
We also identified a number of proteins involved in the metabolism of proteins, carbohydrates, and fats. These metabolisms are related to the acid-base balance in the organism. Acidification can change the acidbase balance in mollusks, and acid-base regulation is closely related to the internal immune defense function.

Effects of OA on proteins associated with biological calcification and mineralization
In this study, for the differentially identified proteins involved in the calcification biological process, in comparison with E1 and C1, troponin T decreased 3.02folds in E1, and the abundance of troponin T protein was decreased in gills of pacific oyster Crassostrea gigas [34]. CaM kinase II alpha decreased 11.59-fold in E2. The biological process of detection of calcium ion was downregulated in E1. In the subcellular localization of identification protein, a portion of the upregulated proteins in E1 was found to be in the calcineurin and troponin complexes in E2; these complexes are closely related to the regulation of calcium.
In previous studies, OA can negatively affect calcification [58]. Larvae of mollusks may be particularly at risk of acidification in marine environments [59] because larval shells are fragile [60]. In OA, biomineralization of gastropods and bivalves is especially vulnerable in mollusk embryonic development and larval growth stage. Biocalcification requires the formation of early skeletal structures, which renders larvae sensitive to the reductions in the saturation state of carbonate ions [61] In calcifying shellfish species, the calcification rate was generally reduced in decreased pH [59]. The mineral formation and calcification mechanism of mollusks and sea urchin larvae also indicate their sensitivity to OA [62]. The decrease of 0.3 pH units resulted in slow growth and calcification of C. gigas during the first 3 days of larval development [63]; the growth and metamorphosis of the C. gigas larvae were significantly reduced in decreased pH [64]. Miller et al. [65] found that the larval shell area of Crassostrea virginica decreased 16% in the range of pH 7.7-8.1; the proportion of calcium in the shells decreased nearly half. OA made larvae shells of the M. galloprovincialis and P. martensiito grow slowly and become deformed [14,48,66]. Ruditapes philippinarum had slightly dissolved calcium shell under the conditions of acute CO 2 acidification [67].
In addition, Dineshram et al. [30] showed that the expression rates of over 24 proteins involved in calcification were changed under decreased pH by comparative proteome analysis, the majority of those proteins were downregulated, including calmodulin and several calcium binding proteins. In previous studies, proteins involved in calcification and acid-base regulation were downregulated under decreased pH in sea urchins [68,69], bivalves [27,28], and other marine invertebrates [70][71][72][73]. The results suggested that acidification can cause the response of calcium-related proteins in the development of the B. areolata larva.

Effects of OA on proteins associated with metabolism
In the present study, the downregulated proteins included proteins, which are involved in carbohydrate metabolic process and ATP binding process in E1. The downregulated proteins in E3 included proteins, which are involved in tricarboxylic acid cycle. In this study, differential proteins involved in energy metabolism included ATP synthase, malate dehydrogenase precursor, arginine kinase, and triosephosphate isomerase.
The result showed that proteins related to the metabolism of the snail larvae have changed under the condition of acidification.
Lower invertebrate animals are easily affected by OA due to their restricted metabolism and physiology [54]. Michaelidis et al. [56] found that the protein synthesis function and respiratory energy of the Mediterranean mussel (M. galloprovincialis) can be reduced if they live in sea water with low pH for a long time. In the pH range of 7.20-8.10, the metabolic substrates are mainly carbohydrates in the thick shell mussel, and no significant difference in the respiration entropy was found; this finding shows that no significant effect on the energy metabolism of the thick shell mussel in this pH range exists. When pH was 7.0, a significant change in energy metabolism was observed. Lannig et al. [74] found that the metabolic pathways of Pacific oysters (C. gigas) have changed under the influence of OA. As marine acidification intensifies, the metabolic pathways of C. farreri also changed [47]. Some species have a strong ability to maintain the acidbase balance. In the face of OA, maintaining the stability of pH value requires a lot of energy; thus, species need to strengthen their metabolism to provide energy. Some species have a weak ability to maintain the acid-base balance. Given the insufficient ability to compensate the effects of pH change, before the appropriate environmental conditions recover, species have to reduce material and energy consumption by reducing the metabolic intensity to maintain their survival [75].
In this study, the proteins involved in protein glycosylation, protein dephosphorylation, lipid metabolism, and ATP synthesis were upregulated in E2. Dineshram et al. [27] found that vacuolar H [+]-ATPase was seen upregulated that indicated that oyster larvae respond to OA by altering its cellular proton transportation at the expense of ATP. In the present study, when E2 was compared with C2, arginine kinase decreased 11.38-fold in E2, arginine kinase catalyzes the reversible transfer from ATP to arginine to form phosphoarginine, there is a negative correlation between arginine kinase abundance and ATP concentration, implying the reduction of energy consumption [34]. In this work, downregulation of malate dehydrogenase (MDH) in E2 (Table S2) perhaps illustrates the influence of elevated pCO 2 exposure on energy metabolism [34]. The proteins involved in protein biosynthesis were upregulated in E3. In this study, acidification pH was about 7.6; the substrate of energy metabolism in B. areolata probably had changed from carbohydrates into protein and lipids with increasing larval development and acidification time; the results need to be further studied and discussed. If the metabolic substrate is converted to protein, then this conversion can lead to break more muscle tissues down for metabolism; this result was detrimental to the healthy growth of mollusks. In the present study, the differentially expressed proteins also included actin, tubulin, and myosin; the result can be related to muscle metabolism.
Following exposure to CO 2 stress, marine calcifiers often exhibit a decrease in metabolic activity in an attempt to conserve energy, which in turn leads to lower calcification rates and somatic growth [36,54,74,76]. Dineshram et al. [30] showed that larvae showed to undergo metabolic suppression in response to decreased pH, while larvae maintained their important cellular functions and metabolic pathways suppressed, including amino acid, carbohydrate metabolism, and deactivation of neuronal signaling. Ko et al. [64] have also shown similar result in larval growth and metamorphosis. These proteins changes reflect a changed energy allocation strategy that may make the B. areolata larvae to survive in response to decreased pH. This study has support from a recent larval proteomics study in the C. gigas) [30].

Other proteins
In this work, proteasome subunit N3 and cathepsin L in E2 (Table S2) were upregulated in E2; ubiquitinconjugating enzyme was upregulated in E1 and E2. Proteasome is mainly responsible for the degradation of abnormal proteins and many cellular regulatory proteins through ubiquitin pathways. Cathepsin L played a major role in protein degradation [77]; in this study, cathepsin L was upregulated in E2 under elevated pCO 2 exposure. The results perhaps indicated that the elevated pCO 2 exposure increased the degradation processes of abnormal proteins.

Conclusions
In the present study, we investigated the proteomic and phenotypic responses of B. areolata larvae in response to OA. This finding suggests that acidification seriously affects veliger larvae after attachment. The identification of differential abundance of key proteins, such as glutaredoxin, heat-shock protein 70, thioredoxin, catalase, cytochrome-c-oxidase, peroxiredoxin 6, troponin T, CaM kinase II alpha, proteasome subunit N3, and cathepsin L during early development, is important to understand molecular mechanisms in the larval developmental of B. areolata in environments with reduced pH. However, further proteomic studies are required to validate these results. Moreover, information about protein abundance must be integrated with preceding patterns of gene expression, such as RT-PCR, and measures of physiological quality to better understand the mechanisms that can permit mollusk larvae to cope with future pH change. Proteome analysis and the novel proteins identified for future studies can extend our understanding of the physiological adaptation of shellfishes to OA.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Fig. S1. pCO 2 manipulation system (Guo et al., 2015). pCO 2 manipulation system with parts as follows 1: air compressor, 2: filter, 3: soda lime column, 4: CaCl 2 anhydrate column, 5: disc-type air filter, 6 and 10: pressure regulation valve, 7 and 11: needle valve, 8 and 12: mass flow sensor, 9: CO 2 cylinder, 13: plastic jar, 14: CO 2 detector, 15: computer. Solid line: air or CO 2 flow; dashed line: digital signal transferring to computer.   Table S1. Differentially expressed proteins with important physiological functions between C1 and E1. Table S2. Differentially expressed proteins with important physiological functions between C2 and E2. Table S3. Differentially expressed proteins with important physiological functions between C3 and E3. Table S4. Differentially expressed proteins with function and process classifications between C2 and E2.