Potential trade-offs between biomineralization and immunity revealed by shell properties and gene expression profiles of two closely related Crassostrea species

Estuaries and coastal zones are areas of high productivity characterized by a high degree of abiotic stress. Many oyster species, including the Pacific oyster (Crassostrea gigas) and the Eastern oyster (Crassostrea virginica), are ecosystem engineers in estuarine and coastal habitats and economically important species supporting multimillion dollar fisheries and aquaculture worldwide (FAO, 2018). Like many benthic sessile species, oysters have a limited ability to escape environmental stressors or predators and rely on mechanical protection (provided by their shells) (Kocot et al., 2016) and physiological and biochemical adjustments (Sokolova et al., 2011) to survive predator attacks and/or shifts in the environmental conditions. Furthermore, because of their filter-feeding habit, oysters are exposed to a large amount and diversity of parasites and pathogens, with an average oyster encountering tens to hundreds of thousands of microbes per second (Allam and Pales Espinosa, 2016; Ben-Horin et al., 2015). In past decades, oyster diseases have expanded in range and increased in severity, often causing staggering losses in the impacted species (Barbosa Solomieu et al., 2015; Beck et al., 2011; Ford and Smolowitz, 2007; Ford, 1996). These conditions put high selective pressure on the protective mechanisms against abiotic and biotic stressors for oyster survival. While oysters are ecologically dominant species in many estuaries, there is considerable variation in the tolerance to disease and environmental stressors between closely related species of oysters. Thus, the Pacific (C. gigas) and Eastern (C. virginica) oysters have similar ecology and feeding habits, yet display a stark difference in disease resistance (Barbosa Solomieu et al., 2015; Elston, 1993; Salvi et al., 2014). Crassostrea gigas is disease tolerant and rarely suffers from epizootic diseases despite serving as a host to a variety of bacteria, viruses and eukaryotic parasites (Barbosa Solomieu et al., 2015; Chu et al., 1996; Elston, 1993). In contrast, C. virginica commonly experiences mass mortalities due to epizootic diseases caused by unicellular eukaryotic parasites Haplosporidium nelsoni and Perkinsus marinus as well as bacterial and viral agents (Barbosa Solomieu et al., 2015; Elston, 1993; Ford, 1996). These remarkable species-specific differences in disease susceptibility between C. gigas and C. virginica have been linked to differences in immunocompetence (Guo et al., 2015). Thus, C. gigas immune cells can better avoid the subversion of the immune defense by the parasite and are more capable of parasite destruction compared with C. virginica hemocytes (Foster et al., 2011; Goedken et al., 2005; Hughes et al., 2010; Sunila and LaBanca, 2003). The hemolymph of C. gigas also contains higher levels of protease inhibitors that prevent cell lysis by the eukaryotic parasite P. marinus compared with the hemolymph of C. virginica (MacIntyre et al., 2003; Romestand et al., 2002). These data show that the interspecific differences in cell-based and humoral immunity may contribute to the differences in disease susceptibility between the Pacific and Eastern oysters. However, the molecular mechanisms underlying these differences are not yet fully understood.

Hemocytes are multifunctional cells that play a key role in innate immune defense, wound repair and biomineralization of oysters (Wang et al., 2018). Oyster hemocytes are involved in recognition, opsonization and destruction of pathogens through phagocytosis, oxidative burst and production of a wide range of antimicrobial compounds (Arumugan et al., 2000; Defer et al., 2013; Gourdon et al., 2001; Hellio et al., 2007). They recognize and encapsulate foreign bodies and participate in wound repair (Ford, 1996; Perez and Fontanetti, 2010). Recent studies have also demonstrated a key role of hemocytes in biomineralization as a result of their involvement in mineral transport as well as in the production of the extracellular matrix (ECM) of the shells (Ivanina et al., 2017; Johnstone et al., 2015; Li et al., 2016; Mount et al., 2004). Notably, different populations of hemocytes of C. gigas (separated by their different floating density) are mostly specialized on either immune defense or biomineralization (Ivanina et al., 2017; Wang et al., 2017b). The developmental mechanisms and physiological consequences of the functional diversity of oyster hemocytes are not yet fully understood, and it remains unclear whether different types of hemocytes represent true differentiated subpopulations, different developmental stages or functional plasticity (Rebelo et al., 2013; Wang et al., 2017b). Regardless of the mechanisms of the functional differentiation, specialization on either immune defense or biomineralization might imply a potential trade-off between these two functions of hemocytes.

The aim of the present study was to investigate the potential interspecific differences in hemocyte involvement in immune defense and biomineralization in two closely related oyster species, C. gigas and C. virginica. We hypothesized that the species-specific differences in the shell properties and disease resistance of oysters are reflected in the prevailing specialization of their hemocytes on either biomineralization or immune defense. To test this hypothesis, we assessed expression of the key biomineralization genes (including genes involved in ion regulation, Ca²⁺ transport, cell–cell interactions and ECM formation) in two main biomineralizing tissues (the hemocytes and the mantle) and assessed the microstructure and mechanical properties of the shells of C. gigas and C. virginica. We also determined mRNA expression of immune defense genes (including molecular pattern-recognition genes and humoral and inflammation-related factors) and the capacity for phagocytosis and mineral transport in the hemocytes of these two congeners. This study provides important insights into the relative immunocompetence and capacity for biomineralization of hemocytes of two ecologically and economically important Crassostrea species that may have implications for mechanistic understanding of the species-specific differences in shell properties, mechanical protection and ability to deal with pathogens.

Oysters C. gigas (Thunberg 1793) (from Fanny Bay, BC, Canada) and C. virginica Gmelin 1791 (from Delaware Bay, NJ, USA) were purchased from a local supplier (Inland Seafood, Charlotte, NC, USA). Oysters were kept in tanks with artificial seawater (ASW; Instant Ocean^®, Kent Marine, Acworth, GA, USA) at 10±1°C and salinity 30±1. They were fed ad libitum with a commercial algal blend containing Nannochloropsis oculata, Phaeodactylum tricornutum and Chlorella spp. (DT's Live Marine Phytoplankton, Sycamore, IL, USA). Algal blend (2–3 ml per 20–25 animals) was added to the experimental tanks every other day.

Hemolymph was extracted from the adductor muscle of oysters using a sterile 10 ml syringe with a 21 gauge needle containing 1 ml of ice-cold salinity 30 ASW [for functional hemocyte analyses, gene expression studies and nanoscale secondary ion mass spectrometry (NanoSIMS) analyses of hemocytes] or 1 ml of ice-cold Alsever's solution containing 20.8 g l⁻¹ glucose, 8 g l⁻¹ sodium citrate, 3.36 g l⁻¹ EDTA, 22.3 g l⁻¹ NaCl (for mineralogical analyses) to prevent aggregation of the hemocytes. For hemocyte functional analyses and mRNA expression, hemolymph was collected from individual oysters. Soft tissues (the gills, the central part of the mantle and the outer mantle edge) were collected from the same oysters. For Fourier-transform infrared spectroscopy (FTIR) spectroscopy and inductively coupled plasma mass spectroscopy (ICP-MS), hemolymph from 7–10 oysters was pooled to obtain sufficient amount of hemocytes. Hemocytes were enumerated using a Brightline hemocytometer and collected by centrifugation for 10 min at 1000 g. For functional studies (phagocytosis and adhesion capacity), hemocytes were analyzed immediately; for all other analyses, hemocytes and soft tissue samples were snap-frozen in liquid nitrogen and stored at −80°C until further analysis.

For the analysis of Ca and Mg content, subpopulations of hemocytes were separated by floating density on a discontinuous Percoll density gradient as described elsewhere (Ivanina et al., 2017). Briefly, isolated hemocytes from 7–15 oysters were resuspended in 4 ml of ice-cold Alsever's solution, layered on a discontinuous Percoll gradient (9.2%, 24.8%, 41.0%, 57.8% v/v in Alsever's solution) and centrifuged at 550 g and +4°C for 40 min. Hemocytes concentrated at the interfaces of the adjacent gradient steps and at the bottom of the tube were collected and labeled H1, H2, H3 and H4 (from top to bottom), respectively. The suspensions containing separate hemocyte subpopulations were diluted 10- to 12-fold with Alsever's solution and centrifuged at 1000 g for 10 min to eliminate Percoll and collect the cells. The pellets were resuspended in Alsever's solution, snap-frozen in liquid nitrogen and stored at −80°C.

The phagocytosis assay was performed as described elsewhere (Ivanina et al., 2014). Briefly, a suspension of Neutral Red-stained, heat-stabilized zymosan (Sigma-Aldrich, St Louis, MO, USA) was added to hemocytes at a final concentration of 200 zymosan particles per hemocyte. Hemocytes were incubated for 30 min at room temperature, centrifuged at 1000 g for 10 min and washed with ASW (salinity 30) to remove extracellular zymosan. Cell-free suspensions of zymosan were used as calibration standards. Hemocytes or zymosan standards were incubated for 5 min with 1% acetic acid in 50% ethanol to extract the Neutral Red and quantified at 550 nm on a microplate spectrophotometer (Multiscan GO, Thermo Scientific, Waltham, MA, USA). Results are expressed as the number of ingested zymosan particles per 100 cells.

Isolated hemocytes (5×10⁵ cells) were placed in 1 ml of ASW (salinity 30) in the wells of a 12-well plate (Costar, Corning, NY, USA) and incubated for 2 h at room temperature as described elsewhere (Ivanina et al., 2014). After incubation, the ASW was collected, the wells were surface-washed with 1 ml of ASW, and the ASW from the incubation and the wash was pooled for each well. The samples were centrifuged for 10 min at 1000 g to collect non-adhered cells. The cells were counted using a Brightline hemocytometer, and the adhesion capacity was expressed as the percentage of adhered hemocytes in each sample.

Total RNA was extracted from hemocytes using a ZR RNA MiniPrep™ kit (Zymo Research, Irvine, CA, USA) according to the manufacturer's instructions. Total RNA from the central part of the mantle, outer mantle edge and gills was isolated using TRI Reagent (Sigma-Aldrich) according to the manufacturer's instructions. RNA concentration was measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Pittsburgh, PA, USA). Single-stranded cDNA was obtained from 0.2 µg of the total RNA using 50 U µl⁻¹ SMARTScribe™ Reverse Transcriptase (Clontech, Mountain View, CA, USA) and 20 µmol l⁻¹ of oligo(dT)₁₈ primers. Transcript expression of target genes was determined by qPCR using a 7500 Fast Real-Time PCR System (Applied BioSystems/Life Technologies, Carlsbad, CA, USA) and SYBR^® Green PCR kit (Life Technologies, Bedford, MA, USA) using gene-specific primers (Table S1). The qPCR reaction mixture consisted of 7.5 µl of 2× SYBR^® Green master mix, 0.3 µmol l⁻¹ each of the forward and reverse primers, 1.5 µl of 10× diluted cDNA template and water to adjust the volume to 15 µl. The reaction mixture was subjected to the following cycling conditions: 10 min at 95°C to denature DNA and activate Taq polymerase and 40 cycles of 15 s at 94°C, 30 s at 60°C and 30 s at 72°C. Serial dilutions of a cDNA standard were amplified in each run to determine apparent amplification efficiency and run-to-run amplification variability (Pfaffl, 2001). The target gene mRNA expression was standardized relative to β-actin mRNA expression (Sanni et al., 2008).

Total hemocyte population and density-segregated hemocyte subpopulations H1, H2, H3 and H4 were used in these studies. Prior to processing, the cell density in each sample was determined using the cell counter slide. The cells were centrifuged (550 g, 40 min, 4°C) and digested in 35% HNO₃ for 24 h at 65°C. The samples were diluted with distilled de-ionized water to achieve the final acid concentration of 3.5% HNO₃ and filtered through 0.20 μm Minisart Plus membrane filters (Sigma-Aldrich). Calcium concentration was analyzed by ICP-MS using a Thermo iCAP 6000 Series (ThermoFisher, Cambridge, UK) and expressed in nmol per cell. Five measurements obtained from the same sample were averaged and used as one biological replicate and three replicates were used for each hemocyte fraction.

Hemocytes from two oysters were pooled, collected by centrifugation and resuspended in 1 ml ASW (salinity 30). Hemocyte suspension (100 μl) was plated on indium–tin–oxide-coated glass slides (Praezisions Glas & Optik GmbH, Iserlohn, Germany) and incubated for 1 h at room temperature in a humid chamber to permit hemocyte attachment. After 1 h, the slide was briefly washed in ultrapure water to remove salts, fixed for 10 min in 70% ethanol, dehydrated by sequential incubation in 95% and 100% ethanol and air-dried for NanoSIMS analysis. For each species, three biological replicates (each consisting of the pooled hemocytes of two individuals) were analyzed. Two additional biological replicates of C. gigas were found to be on slides without sufficient conductibility. They were therefore sputter-coated with gold (Cressington 108 auto sputter coater) and reanalyzed. SIMS imaging was performed using a NanoSIMS 50L instrument (Cameca, Gennevilliers, France). Three individual frames were analyzed for each biological replicate.

The Cs source of the NanoSIMS instrument was used to erode and ionize atoms of the sample. Among the received secondary ions, images of ¹²C⁻, ¹⁶O⁻, ¹²C¹⁴N⁻ and secondary electrons were recorded simultaneously for areas with sufficient cell density. Prior to the analysis, sample areas of 50×50 µm were sputtered for 90 s at 600 pA to reach the steady state of secondary ion formation and clean the surface. The primary ion beam current during the analysis was 1 pA; the scanning parameters were 512×512 pixels for areas of 47×47 µm, with a dwell time of 250 µs per pixel; 60 planes were analyzed. Slightly smaller areas (28–35 µm) of the same locations were then analyzed with the NanoSIMS O source (Duoplasmatron) and ²³Na⁺, ²⁴Mg⁺ and ⁴⁰Ca⁺ ions were recorded. It is worth noting that SIMS analyses are destructive; therefore, the sequential analyses with the Cs and the O source deliver signals from different sample depths, which might bias the Cs/O correlations. However, the sample erosion by the Cs source was shown to be small enough to resolve cell size features (5.7 nm for a similar primary Cs beam in a smaller area; Saka et al., 2014). The primary ion beam current during the analysis was 25 pA; the scanning parameters were 512×512 pixels with a dwell time of 250 µs per pixel; 60 planes were analyzed. Secondary ions were detected using mass detectors equipped with electron multipliers (Hamamatsu, Hamamatsu City, Japan). The instrument was adjusted to reach a mass-resolving power appropriate to suppress interference at all masses, allowing the separation of ¹²C¹⁴N⁻ from interfering ions such as ¹²C₂¹H₂⁻.

NanoSIMS data were analyzed using Look@NanoSIMS software (Polerecky et al., 2012). Depending on the signal quality, 58–60 planes were drift corrected and accumulated. The cells were defined as regions of interest (ROIs) employing the Na⁺ signal and verified by comparison with the CN⁻ and secondary electron pictures. For the gold-coated samples, the accumulated O⁻ and Na⁺ were used to align the two measurements. To avoid potential bias due to differences in the analytical procedure, the gold-coated samples were excluded from the quantitative analyses of ion ratios. For all other samples, Ca/Mg, Ca/Na and Mg/Na ratios were calculated for each ROI representing individual cells. Cell sizes were calculated from the number of pixels in the ROI assuming an elliptical cell shape. Ca/Na and Mg/Na ratios were used as proxies for the relative Ca and Mg concentrations in the hemocytes, assuming that Na density per unit cell area is similar in the osmoconforming mollusks maintained under the same salinity conditions.

Thermal characteristics of samples were measured with a TGA Q500 thermal analysis system (TA Instruments, New Castle, DE, USA). Approximately 10 mg of powdered dry oyster shell was used for each test. All TGA experiments were performed under nitrogen atmosphere (balance gas: N₂ purge rate 40.0 ml min⁻¹; sample gas: N₂ purge rate 60.0 ml min⁻¹), from 25 to 900°C, at a heating rate of 10°C min⁻¹. Three samples per species were analyzed.

Blocks of 1.0×3.0 cm were cut from the central portions of the right valves and weighed using a Mettler Toledo XP26DR Microbalance (Mettler Toledo, LLC, Columbus, OH, USA). The volume of the samples was assessed using the liquid displacement technique in 5 ml graduated cylinders. Five samples per species were measured.

Hemocytes isolated from the pooled hemolymph of 10 oysters were resuspended in 4 ml Alsever's solution and layered onto 10 ml of 10% sucrose in Alsever's solution. Cells were purified by centrifugation for 30 min at 550 g, 4°C. The resulting pellet was resuspended in 1 ml of 10% sucrose–Alsever's solution, vigorously vortexed for 20 s, mixed with 4.0% NaOCl (1:2 v/v hemocyte suspension to NaOCl ratio), and centrifuged for 90 s at 7000 g. The supernatant was removed and replaced with 2.5% NaOCl. After 20 min incubation at room temperature, the samples were centrifuged for 5 min at 8000 g and rinsed with ultrapure water adjusted to pH 10 with NaOH. The tubes were flash frozen in liquid nitrogen and freeze-dried overnight (Labconco Freezone 4.5, Kansas City, MO, USA). Dry samples were mixed with KBr, pressed into pellets and studied using a FTIR spectrometer (Broker Vertex 70, Billerica, MA, USA) in transmission mode. For shells, material from the central portions of the left valves of five oysters of each species were ground and pressed into KBr pellets. Pellets were analyzed using a FTIR spectrometer (Bruker Vertex 70) in transmittance mode with a resolution of 4 cm⁻¹ and 128 scans per sample. Absorbance spectra were processed using Spectrum 5.1 software (Perkin Elmer, Shelton, CN, USA). The 600–1000 cm⁻¹ regions of the spectra were isolated, baseline corrected and normalized, and the positions and heights of absorption peaks of ν4 and ν2 vibrational bands were measured. The ν4/ν2 band intensity ratio is inversely proportional to calcite crystal and is commonly used as a proxy for the crystallinity of the mineral phase (Beniash et al., 1997; Gueta et al., 2007).

The samples were prepared in a similar way to the procedure for FTIR spectroscopy. Dry samples were resuspended in a small volume of ethanol and mounted on carbon-coated Ni grids. The samples were studied using a Jeol 2200 (Akishima, Tokyo, Japan) high-resolution transmission electron microscope, equipped with a field emission gun, in bright-field and SAED modes at an accelerating voltage of 200 kV.

Blocks (1 cm×1 cm) were cut from the central portion of the right valves from five C. gigas and five C. virginica oysters with a shell length of 6–8 cm. The blocks were mounted in Epofix fast-curing epoxy resin (EMS, Hatfield, PA, USA) and left to polymerize for 24 h at room temperature. The blocks were trimmed with a water-cooled slow-speed saw (IsoMet, Buehler, Lake Bluff, IL, USA) and polished in the transverse plane along the long axis of the shell using Metadi diamond suspensions down to 0.25 µm (Buehler) in a saturated CaCO₃ solution. The microhardness tests were conducted using a microindenter (IndentaMet 1104, Buehler) at 0.245 N load and 5 s dwelling time using Vickers indenter tip. Indentations were made in the middle of the shells' foliated layer. Five indentations were made per sample at least 100 µm away from other indents. Data from five indentations were averaged for each shell sample to obtain a biological replicate.

Blocks of shell mineral (1.5 cm×4.5 cm) were cut from the central portion of the right valves from five C. gigas and five C. virginica oysters with a total shell length of 6–8 cm. The three-point bending tests were conducted with an electronic universal mechanical testing machine (Instron 5566, Norwood, MA, USA). The applied velocity of the bending load was 1 mm min⁻¹, with 5 kN static load, 0.5 N pre-load and 22 mm support span. Five measurements of stiffness obtained from the same shell were averaged and used as a biological replicate.

Blocks of shell (1 cm×1 cm) were cut from the central portion of the right valves from five C. gigas and five C. virginica oysters, mounted in Epofix fast-curing epoxy resin (EMS, Hatfield, PA, USA) and allowed to polymerize for 24 h at room temperature. The blocks were trimmed with a water-cooled slow-speed saw (IsoMet, Buehler), polished in the transverse plane along the long axis of the shell using Metadi diamond suspensions (Buehler) down to 0.25 µm diamond particle size in saturated CaCO₃ solution, and sputter coated with gold under vacuum (Anatech Hummer 10.2, Smyrna, DE, USA). The samples were examined by field-emission SEM (Zeiss Sigma 500 VP FE-SEM, Oberkochen, Germany, with Oxford Aztec X-EDS, Oxford, UK) at 15 kV.

One- and two-way ANOVA were used to test the effects of species on the studied traits (Tables 1 and 2). For ion ratios based on the NanoSIMS analyses, cell size was included as a covariate to correct for the potential cell size and/or margin effects. The average cell size of the NanoSIMS samples was similar for the two studied species (P>0.05) (4.2±0.7 µm, N=19 and 4.6±0.4 µm, N=52 for C. gigas and C. virginica, respectively). Prior to analyses, data were tested for normality and homogeneity of variance by Kolmogorov–Smirnoff and Levene's tests, respectively, and normalized as needed using the Box–Cox common transforming method. Fisher's least significant differences (LSD) tests were used for planned post hoc comparisons of the differences between pairs of means of interest. Differences were considered significant if the probability of Type I error was less than 0.05. The data are presented as means±s.e.m. unless indicated otherwise. Pearson correlation analysis was conducted to assess correlations between mRNA levels of target genes, tissue types and species. Principal component analysis (PCA) was used to reduce the dimensionality of the dataset and determine the groups of molecular traits (gene transcription profiles) that distinguish the different tissue types and/or species. The correlation analysis and the PCA were conducted on raw (non-transformed) data. All statistical analyses were performed with the OriginPro 2015 software package (Northampton, MA, USA).

Tissue-specific expression of carbonic anhydrase (CA) isoforms differed between the two studied congeners (Fig. 1A–C, Table 1). CAI was expressed at relatively low levels (0.5–0.6 CAI to β-actin ratio) in all studied tissues of C. gigas (Fig. 1A). In C. virginica, CAI mRNA expression was low in hemocytes and the mantle, but considerably higher (∼8.3 CAI to β-actin ratio) in the gill (Fig. 1A). When compared between the two studied species, CAI mRNA levels were considerably higher in the gills of C. virginica compared with those of C. gigas (Fig. 1A). In contrast, CAXIV was highly expressed in the mantle of C. gigas (52–86 CAI to β-actin ratio), exceeding the respective levels in hemocytes and the gill by 40- to 66-fold (Fig. 1C). In C. virginica, CAXIV also tended to be more highly expressed in the mantle than in other tissues (∼3 CAXIV to β-actin ratio in the mantle versus 0.5–1 in hemocytes and the gill) (Fig. 1C). When compared between the two species, CAXIV mRNA levels were significantly higher in the mantle of C. gigas compared with that of C. virginica (Fig. 1C). mRNA expression of CAII was similar in all studied tissues and between C. gigas and C. virginica (Fig. 1B).

V-type H⁺-ATPase mRNA showed similar expression levels in all studied tissues of C. gigas (Fig. 1D). In C. virginica, V-type H⁺-ATPase mRNA expression in hemocytes was 2- to 4-fold higher than in other studied tissues (P<0.05). The expression levels of V-type H⁺-ATPase mRNA in the hemocytes and the gills were higher in C. virginica than in C. gigas. Expression of V-type H⁺-ATPase was similar in the mantle of C. gigas and C. virginica (Fig. 1D).

In C. gigas, mRNA levels of Ca²⁺-ATPase Type 2C were significantly elevated in the central mantle compared with other tissues, while in C. virginica the mRNA levels of this gene did not differ in different tissues (Fig. 1E). Expression of Ca²⁺-ATPase Type 2C mRNA was higher in the gill (but not in other tissues) of C. virginica than of C. gigas (Fig. 1E). In contrast, mRNA levels of the plasma membrane Ca²⁺-ATPase were consistently higher in C. virginica than in C. gigas in all studied tissues (Fig. 1F).

In C. virginica, Na⁺/H⁺ antiporter 9 (NHX9) and Na⁺/H⁺ antiporter 3 (NHE3) mRNA expression in the gill was higher than that in other tissues of this species (albeit this trend was only significant for NHE3, P<0.05), whereas in C. gigas the mRNA levels of these two genes showed little tissue-to-tissue variation (Fig. 1G,H). NHX9 and NHE3 transporters tended to have higher expression levels in C. virginica tissues than in those of C. gigas (Fig. 1G,H). This species-specific difference was significant in the gill, where NHE9 and NHE3 expression was ∼27 and ∼7 times higher in C. virginica than in C. gigas, respectively. Overall, the interspecific differences in mRNA expression of genes involved in the ion and acid–base regulation indicate higher activity of the membrane Ca²⁺, Na⁺ and H⁺ transporters in C. virginica than in C. gigas.

In C. gigas and C. virginica, mRNA expression of silk-like protein (SLP) was highest at the mantle edge (Fig. 2A). When compared between the species, mRNA levels of SLP and nacrein were notably higher in the mantle edge of C. virginica versus C. gigas (Fig. 2A,B). Across all tissues, the expression of fibronectin Prot3L was lower in C. gigas than in C. virginica, while the expression of fibronectin Prot2L was lower in C. virginica (Fig. 2C,D). Fibronectin ankyrin showed higher levels of mRNA expression in the mantle and/or the gill versus the hemocytes in both studied species of oyster. Compared between the species, fibronectin ankyrin expression was significantly higher in all tissues of C. virginica (except the mantle edge, where this trend was not statistically significant) (Fig. 2E). Overall, mRNA levels of SLP, fibronectin Prot3L, fibronectin ankyrin and nacrein were higher in the biomineralizing tissues (the mantle and/or hemocytes) of C. virginica than in those of C. gigas. In contrast, the mRNA expression of fibronectin Prot2L was higher in the mantle of C. gigas versus C. virginica (Fig. 2, Table 1).

When compared among different tissues, VEGF and VEGF receptor (VEGFR) mRNA levels were higher in the mantle and the gills than in hemocytes of C. virginica; a similar, albeit less pronounced, trend was found in C. gigas (Fig. 3A,B). Overall, mRNA expression of VEGF and VEGFR was higher in the mantle and the gills of C. virginica than in the respective tissues of C. gigas (Fig. 3A,B).

Expression levels of casein kinases (CKI and CKII) involved in ECM protein phosphorylation (Rusenko et al., 1991; Sfeir and Veis, 1996; Veis et al., 1997) were relatively uniform across different tissues within the same species (except for elevated levels of CKI in the gill of C. virginica) and were consistently higher in C. virginica than in C. gigas (Fig. 3). In contrast, mRNA levels of chitin synthases (CSII and CSIII) were higher in the biomineralizing tissues (the mantle and/or the hemocytes, respectively) (Fig. 3). When compared between the species, CSII mRNA levels were higher in the mantle edge of C. gigas versus C. virginica. Expression of CSIII was higher in the hemocytes and lower in the mantle in C. gigas versus C. virginica (Fig. 3).

Hemocyte concentrations in oyster hemolymph were slightly but significantly higher in C. gigas (11.9±0.7 million cells per ml, N=10) than in C. virginica (9.7±0.6 million cells per ml, N=10) (P<0.05). Phagocytic activity was similar in C. gigas and C. virginica hemocytes (980±11 versus 820±8 zymosan particles per 1000 cells, N=10 and 9, respectively) (P>0.05). The adhesion capacity of hemocytes was high in both studied species with >96% hemocytes attached after 2 h of incubation.

Of the six studied molecular pattern-recognition genes, five (TLR2, TLR3, TLR4, mannose receptor and killer cell lectin-like receptor) tended to be more highly expressed in C. gigas than in C. virginica (Fig. 4, Table 2). Among these, mRNA levels of TLR3 and killer cell lectin-like receptor were significantly higher in hemocytes of C. gigas (by ∼7- and ∼45-fold, respectively) than in C. virginica (P<0.05) (Fig. 4B,E, Table 2). mRNA expression of TLR2, TLR4, mannose receptor 2 and scavenger receptor cytosine rich (SRCR) did not significantly differ between hemocytes of the two studied species (Fig. 4A,C,D,F).

Expression levels of mRNA encoding big defensin, complement component 3 (Cq3) and tumor necrosis factor (TNF) were ∼433-, ∼6- and ∼4-fold higher, respectively, in the hemocytes of C. gigas versus C. virginica (Fig. 4H–J). The mRNA expression of C-type lectin was similar in the hemocytes of the two studied species (Fig. 4G). The mRNA expression level of integrin was similar in the two studied species (P=0.07, Table 2).

PCA of species- and tissue-specific variation in the expression of biomineralization-related genes showed that the species had a relatively strong loading in the first principal component (PC1), which explains 41% of variation (Fig. 5A; Table S2). The second principal component (PC2), explaining 15% of variation, had a very strong loading associated with tissue type. Interestingly, the two mantle locations (mantle edge and central mantle) clustered together, while the gills clustered with the hemocytes (Fig. 5A; Table S2). Pearson correlation analysis revealed that higher gene expression levels of multiple genes was associated with C. virginica (Table S3). Specifically, statistically significant positive correlations were found between expression of 10 out of 18 studied biomineralization-related genes in C. virginica, while higher expression of only two genes (CA14 and fibronectin Prot2L) was associated with C. gigas.

When only hemocytes were considered in the analysis, species showed a strong loading in PC1, which explains 44% of variation (Fig. 5B; Table S4). PCA analysis of the expression of biomineralization-related and immune genes showed a clear separation of C. gigas and C. virginica hemocytes, with higher expression of the molecular pattern-recognition genes (except SRCR) as well as humoral and inflammation-related factors in hemocytes of C. gigas versus C. virginica (Fig. 5B; Table S3). High expression of cell matrix-interacting genes such as β-integrin and fibronectin ankyrin was associated with C. virginica. Moreover, expression of CKs and ionoregulatory genes was positively associated with C. virginica (Fig. 5B). Based on the Pearson correlation analysis, expression of five genes tentatively associated with biomineralization, pH control and ion transport (CK1, VEGFR, NHX9, PM Ca²⁺-ATPase, V-type H⁺-ATPase) was strongly and positively correlated with C. virginica (Table S5). Expression of the genes potentially associated with defense (TLR3, Cq3, TNF, defensin) and cell adhesion (C-type lectin and fibronectin Prot2L) positively correlated with C. gigas (Table S5).

Calcium content was significantly lower in the total hemocyte population of C. gigas than in that of C. virginica (0.41±0.028 versus 0.72±0.09 nmol cell⁻¹; P=0.02) (Fig. 6A). This was also reflected in the significantly lower Ca content of individual fractions (H1–H4) of C. gigas hemocytes compared with values for C. virginica. In C. gigas as well as C. virginica, Ca content of hemocytes increased from the lightest to the heaviest hemocyte fraction, and this trend was more pronounced in C. gigas (Fig. 6A).

FTIR analysis revealed that the mineral phase from the hemocytes of C. gigas and C. virginica consists of calcium carbonate, predominantly calcite (Fig. 6B,C). Notably, different preparations of the mineral phase from replicate samples of hemocytes from the same species showed variations in the spectra of the mineral (Fig. S1). Our HR-TEM/SAED studies identified multiple mineral phases of CaCO₃ in hemocytes, including amorphous calcium carbonate (ACC), calcite, dolomite and vaterite (Fig. 7). These data indicate that the calcium carbonate mineral inside oyster hemocytes is highly unstable and potentially amorphous, and that it spontaneously transforms into other calcium carbonate phases (such as calcite) during the preparation procedure.

Analysis of the NanoSIMS maps showed that the ratios of the intracellular Ca/Na and Mg/Na signals (indicative of the cellular Ca and Mg content, respectively) were ∼2.7- to 3-fold lower in hemocytes of C. gigas than in those of C. virginica (P<0.001 and P=0.003 for Ca/Na and Mg/Na ratios, respectively) (Fig. 8). Larger hemocytes contained relatively more Ca and Mg compared with smaller ones, as demonstrated by the significant positive effects of cell size (used as a covariate in ANCOVA) on Ca/Na and Mg/Na ratios (P=0.002 and P<0.001, respectively). The ratio of the Ca/Mg signal was higher in larger hemocytes (P<0.001 for cell size used as a covariate in ANCOVA) but not significantly different between the species (P>0.05). NanoSIMS images of oyster hemocytes with overlapping hotspots of Ca and Mg (Fig. 8) and Ca and O (Fig. S2) are consistent with the presence of CaCO₃ mineral in oyster hemocytes.

Overall, Ca content analysis revealed that hemocytes of C. virginica carry more calcium and that in both species the heaviest (H4) hemocyte fraction contains more Ca than the others. Furthermore, the results of FTIR and diffraction analyses reveal that the mineral fraction of oyster hemocytes contains several mineral phases of CaCO₃.

FTIR analysis of the oyster shells indicated that in both studied species the major mineral fraction of the shell is calcite, as indicated by the positions of ν₄ (712 cm⁻¹), ν₂ (878 cm⁻¹) and ν₃ (1420 cm⁻¹) absorbance bands (Andersen and Brecevic, 1991) (Fig. 9A,B). The ν₂/ν₄ peak ratio (a measure of the crystallinity of calcite; Beniash et al., 1997; Gueta et al., 2007) was significantly higher in C. gigas than in C. virginica shells (3.5±0.30 versus 2.2±0.07; P=0.002) indicating that the mineral in C. gigas shells is less crystalline than in those of C. virginica.

SEM analysis of the foliated layer revealed similarity of the overall ultrastructure of the shells of C. gigas and C. virginica (Fig. S3). However, the thickness of the laths in the shells of C. gigas was significantly lower than that in C. virginica (187±22 versus 251±24 nm; P=0.03) (Fig. 9C). Similarly, the microhardness of the foliated layers of the shells was significantly lower in C. gigas than in C. virginica [152±9 versus 206±1.3 Vickers Pyramid Number (HV); P<0.005] (Fig. 9D). The three-point bending test revealed that the shells of C. virginica were significantly stiffer than those of C. gigas [24(±4.1)×10⁵ versus 4.2(±2.2)×10⁵ N m⁻¹; P<0.0001] (Fig. 9E). Similarly, the maximum compressive load was lower in C. gigas shells than in those of C. virginica [2.5(±0.56)×10⁵ versus 9.35(±5.5)×10⁵ N m⁻¹; P=0.05] (Fig. 9F). Taken together, these data reveal differences in shell mineralogy and ultrastructure between the two species studied and show that the shells of C. gigas are mechanically inferior to those of C. virginica.

TGA analysis did not show any significant differences in the organic content of the shells between C. gigas and C. virginica (1.95±0.08% versus 3.12±1.01%; P=0.18) (Fig. 9G,H). Notably, the density of the C. gigas shells was significantly lower than that of C. virginica (2.0±0.59 versus 3.2±0.37 g ml⁻¹, P=0.005) (Fig. 9I). Taken together, the results of TGA and density analyses suggest that the shells of C. gigas are more porous than those of C. virginica.

Despite similar shell morphology, the mineralogical, compositional and mechanical properties of the shells dramatically differ between the two closely related species of oyster, C. virginica and C. gigas. FTIR analysis revealed that calcite in the shells of C. gigas is less crystalline than in C. virginica suggesting that the kinetics of crystal maturation process in these shells is slower or that they contain more ACC (Beniash et al., 1997; Gueta et al., 2007). These differences in mineral crystallinity can affect the mechanical properties of the crystals and influence shell hardness. This suggestion is borne out by the results of the mechanical analyses, showing considerably lower shell microhardness in C. gigas than in C. virginica. It is worth noting that the shell microhardness in oysters is a plastic trait and can change depending on the growing conditions. Our earlier studies (Beniash et al., 2010) showed that the shell microhardness of C. virginica juveniles grown under reduced pH (elevated CO₂) was significantly lower than that of juveniles grown under normal conditions. Based on the shell microstructural analysis, we hypothesized that the differences in the microhardness of the shells grown at normal and elevated CO₂ levels are associated with differences in lath thickness (Beniash et al., 2010). In this earlier study, the shells with higher microhardness were made of thinner laths (Beniash et al., 2010). In contrast, in the present study, the shells with thicker laths were harder. These results suggest that it is not the lath thickness (as we speculated earlier) but other factors that affect the microhardness of the lath layers. These differences can be due to the interfacial properties of the laths or the adhesive properties of the organic matrix. It is also possible that the greater hardness of C. virginica shells is associated with the higher mineral crystallinity as assessed by FTIR. Bulk mechanical properties of the shells such as stiffness and maximum load were also significantly greater in C. virginica than in C. gigas. We attribute these differences to the higher porosity of C. gigas shells. Taken together, these data demonstrate that the shells of C. virginica are mechanically superior to those of C. gigas on both the microscale and macroscale.

Biomineralization of oyster shells is a complex process involving the coordinated action of multiple organs (Gardner et al., 2011; Ivanina et al., 2017; Johnstone et al., 2015; Li et al., 2016). Our data indicate stronger expression of acid–base and ionoregulatory genes in different tissues (including mantle, gills and hemocytes) of C. virginica compared with C. gigas. This includes several isoforms of carbonic anhydrases (except CAXIV, which is more highly expressed in the mantle of C. gigas), an ATP-dependent proton pump (V-type H⁺-ATPase) and sodium–proton exchangers that play a crucial role in acid–base homeostasis by regulating intracellular and extracellular pH and bicarbonate concentrations (Larsen et al., 2014). Notably, the transcripts of many acid–base and ion regulation-related genes were as (or more) highly expressed in the gill as in the biomineralization-specialized tissues such as the mantle or hemocytes. This is consistent with the important role of the gill in systemic acid–base and ion homeostasis (Larsen et al., 2014). Furthermore, the plasma membrane Ca²⁺-ATPase involved in intracellular Ca²⁺ homeostasis is also more highly expressed in the different tissues of C. virginica compared with C. gigas. Higher expression of the genes involved in acid–base regulation and Ca²⁺ transport may help create a favorable milieu for biomineralization due to a better regulation of pH and CaCO₃ saturation levels in the extrapallial cavity and at tissue–shell interface of C. virginica versus C. gigas (Holcomb et al., 2014; Ramesh et al., 2017; Venn et al., 2013).

Expression of the ECM-related genes shows complex profiles among different tissues in the two studied species. Generally, C. virginica expressed higher mRNA levels of fibronectins (fibronectin ankyrin and fibronectin Prot3L but not Prot2L) as well as genes associated with ECM maturation (CKs) and SLP in the mantle and/or the gill tissues. The tissue-specific patterns of expression of CSs and nacrein were more variable between the two studied tissues. Interestingly, transcripts of some ECM-related proteins (including fibronectins and CKs) were highly abundant in the gills of oysters, an organ not traditionally associated with biomineralization. This finding is consistent with earlier reports of de-localized synthesis of ECM proteins in oysters and their subsequent transport in the blood to the biomineralization sites (Wang et al., 2013). This might indicate that the biomineralization in mollusks has a greater contribution from organ systems other than the mantle, thus challenging the classical biomineralization model.

Notably, the genes involved in cell–cell signaling (VEGF and VEGFR) showed higher transcript levels in mantle and gill tissues of C. virginica versus C. gigas. VEGF and its receptor have been implicated in biomineralization of invertebrates such as sea urchins, where VEGF regulates migration of the skeletogenic primary mesenchymal cells (Adomako-Ankomah and Ettensohn, 2014; Duloquin et al., 2007) and growth of the calcitic spicules (Knapp et al., 2012). In C. gigas, subpopulations of hemocytes specialized on biomineralization also expressed higher VEGF and VEGFR mRNA levels compared with the immune-specialized hemocytes, thus underscoring the potential role of VEGF in biomineralization (Ivanina et al., 2017). However, our recent study in C. gigas juveniles showed that while chemical interference with VEGF signaling modulates biomineralization, VEGF inhibition does not result in disrupted skeletogenesis (Ivanina et al., 2018). While higher expression of VEGF and VEGFR in C. virginica compared with C. gigas is consistent with the pattern found for most studied biomineralization-related genes, the functional implications of VEGF for oyster biomineralization remain speculative because of the poorly defined role of this signaling pathway in the mollusks.

Our ICP-MS and NanoSIMS studies revealed large differences in the Ca content of the different fractions of hemocytes in both species, with the largest and heaviest hemocytes containing the highest amount of Ca. Notably, in our earlier study we found that the free Ca²⁺ concentration was lowest in the heaviest hemocyte fraction (Ivanina et al., 2017). These results suggest that this hemocyte fraction contains bound Ca potentially in the form of CaCO₃ mineral. Importantly, Ca content was significantly higher in the total hemocyte populations, as well as across different hemocyte size fractions, in C. virginica than in C. gigas. These results are in a good agreement with the results of our transcriptome analysis, showing that the expression of the biomineralization-related genes (especially those involved in ion regulation and Ca transport) is higher in hemocytes of C. virginica. This might imply that biomineralization processes in C. virginica are more active than those in C. gigas, thereby contributing to the superior shell quality of C. virginica versus C. gigas.

The ability of oyster hemocytes (especially those of C. virginica) to accumulate significant amounts of Ca without the associated Ca²⁺ toxicity may be achieved by storing intracellular Ca²⁺ in CaCO₃ mineral particles. Our FTIR spectroscopy and HR-TEM analyses confirmed the presence of calcium carbonate mineral particles in hemocytes. Importantly, several CaCO₃ mineral phases were identified in hemocytes, including the least thermodynamically stable ACC, suggesting that intracellular CaCO₃ is stored as highly unstable ACC, which is transformed into other phases during experimental preparation and/or deposition on the growing shell. The presence of CaCO₃ in oyster hemocytes demonstrated in the present study is consistent with the earlier reports of CaCO₃ in bivalve hemocytes assessed by microscopy (Li et al., 2016; Mount et al., 2004). The biomineralization strategy involving mineral ion sequestration by cells away from the mineralization site coupled with the transport of these ions in the form of a metastable mineral phase were identified in vertebrates (Akiva et al., 2015; Kerschnitzki et al., 2016) and echinoderms (Vidavsky et al., 2014, 2015). Our results suggest that this strategy might be more widespread than was previously thought.

Our study showed considerable differences in the gene expression patterns and functional characteristics of hemocytes of the two closely related species of oyster, indicating that the hemocytes of C. gigas are more strongly specialized on immune defense, while C. virginica hemocytes possess more active ionoregulatory and Ca transport mechanisms.

The pathogen-recognition system of C. gigas hemocytes is more highly expressed than that in C. virginica. Immune surveillance in oysters, as in other organisms lacking adaptive immunity, entirely depends on the germline-encoded pattern recognition receptors (PRRs) that recognize pathogen-associated molecular patterns (PAMPs) from various microbes (Ward and Rosenthal, 2014). Crassostrea gigas, the only Ostreidae species for which full genome data are available, shows massive expansion and diversification of the PRR genes compared with vertebrates, insects or cnidarians, likely reflecting adaptation to the abundant and diverse microbial challenges in estuarine and coastal habitats (Zhang et al., 2015). The present study shows that circulating hemocytes of C. gigas express higher levels of several key PRR-encoding genes compared with the closely related but less disease-resistant C. virginica (Barbosa Solomieu et al., 2015; Elston, 1993; Ford, 1996). This includes three isoforms of TLR receptors (TLR2, TLR3 and TLR3) belonging to a family of transmembrane receptors specialized for the recognition of viral and bacterial PAMPs (Ward and Rosenthal, 2014; Zhang et al., 2015). Functional studies showed that oyster TLRs (including TLR2, TLR3 and TLR4) are strongly activated by bacterial PAMPs such as the heat-killed Vibrios, bacterial lipopolysaccharide and lipoproteins (Zhang et al., 2013). The oyster TLR3 also responds to dsRNA (albeit more weakly than to bacterial PAMPs) and thus may also be involved in virus recognition (Zhang et al., 2013). Interestingly, elevated expression of TLRs in C. gigas hemocytes goes hand-in-hand with high expression of TNFα, a major TLR-dependent regulator of inflammation (Parameswaran and Patial, 2010). Elevated expression of TLRs and TNFα in C. gigas may thus indicate higher activity of the pathogen-induced inflammatory pathways in this species compared with C. virginica and may contribute to higher resistance to bacterial pathogens in C. gigas.

Hemocytes of C. gigas are better equipped to fight eukaryotic pathogens than hemocytes of C. virginica. Crassostrea gigas hemocytes present considerably higher mRNA expression of PRRs involved in the recognition of eukaryotic pathogens including mannose receptor 2 (∼1.5-fold higher than in C. virginica) and killer cell lectin-like receptor (∼45-fold higher than in C. virginica). Mannose receptor 2 and killer cell lectin-like receptor belong to the family of C-type lectins with broad glycan specificity recognizing fungal PAMPs as well as other eukaryotic pathogens including protists and parasitic helminths (Llibre et al., 2016; van Die and Cummings, 2017; Vautier et al., 2012). Taken together, the species-specific transcription patterns of immune-related genes suggest that C. gigas hemocytes are better equipped to recognize a broad repertoire of bacterial, viral and eukaryotic pathogens compared with those of C. virginica.

The antibacterial humoral defense of hemocytes also appears to be more active in C. gigas than in C. virginica. Remarkably, the basal mRNA expression levels of a major humoral component of the immune system, big defensin, was ∼432-fold higher in C. gigas hemocytes than in C. virginica. Big defensin is an important antimicrobial peptide (AMP) with a broad pathogen specificity and microbicidal activity against Gram-positive and Gram-negative bacteria and eukaryotic microbes (Rosa et al., 2011; Schmitt et al., 2012). In C. gigas, big defensins are highly inducible by exposure to bacterial infections and/or to PAMPs (Rosa et al., 2011; Schmitt et al., 2012). This inducibility, combined with the high baseline levels of big defensin, may aid C. gigas hemocytes in effective elimination of bacterial and protozoan pathogens. Furthermore, C. gigas hemocytes express considerably (∼6-fold) higher mRNA levels of Cq3 compared with C. virginica hemocytes. Cq3 acts as a central innate immunity hub mediating opsonization and production of AMPs to eliminate pathogens and activating the complement cascade as the first line of immune defense (Ricklin et al., 2016; Volanakis, 1990). In marine bivalves such as the mussel Mytilus corruscus, the razor clam Sinonovacula constricta and oyster C. gigas, Cq3 is stimulated by bacterial infection (Vibrio and Micrococcus species), exposure to bacterial PAMPs or toxic insults (such as dissolved copper) (Chen et al., 2018; Dong et al., 2017; Peng et al., 2017, 2016; Wang et al., 2017a). A highly expressed Cq3 system combined with the elevated expression of AMPs such as big defensin might contribute to a more rapid activation of the humoral immune defense and more effective pathogen killing in C. gigas compared with C. virginica.

Transcriptome analysis also sheds light on the mineral sequestration and transport mechanisms by hemocytes, showing strong species-specific differences in the expression of genes involved in ion regulation and deposition of the ECM proteins. Thus, the transcripts of the genes involved in Ca²⁺, Na⁺ and H⁺ transport (V-type H⁺-ATPase, Ca²⁺-ATPase, plasma membrane Ca²⁺-ATPase and NHE9 sodium–proton exchanger) were consistently more highly expressed (by ∼2- to 5-fold) in the hemocytes of C. virginica versus C. gigas. Active ionoregulation is crucially important for bivalve hemocytes as they sequester and transport Ca²⁺ to the sites of biomineralization (Ivanina et al., 2017; Mount et al., 2004; Wang et al., 2017c). This poses major physiological challenges to hemocytes as they need to carry large quantities of Ca²⁺ to the sites of shell biomineralization, while at the same time avoiding cytotoxicity. This might explain the higher expression of genes involved in the maintenance of intracellular Ca²⁺ levels (notably, plasma membrane Ca²⁺-ATPase) (Brini and Carafoli, 2011). In vertebrates and Drosophila, plasma membrane Ca²⁺-ATPase is mostly responsible for the efflux of Ca²⁺ from the cells to regulate the household Ca²⁺ homeostasis (Brini and Carafoli, 2011). If similar mechanisms exist in oysters, high expression of plasma membrane Ca²⁺-ATPase in C. virginica hemocytes may aid in the regulation of intracellular free Ca²⁺ levels of the cells that accumulate and transport Ca to the biomineralization sites. Transcripts of the genes encoding CKs involved in the phosphorylation of ECM proteins were ∼7–10 times more abundant in C. virginica hemocytes. In contrast, genes encoding ECM proteins (including SLP, nacrein and fibronectins) were similarly expressed in the hemocytes of the two studied species of oyster.

Our study demonstrates the molecular and cellular basis of the species-specific differences in important fitness-related phenotypes (disease resistance and exoskeleton properties) between two closely related and ecologically similar species of oyster. The hemocytes of the two studied species are multifunctional cells, yet transcriptomic and mineralogical profiles of the hemocytes show predominant specialization on immune defense (C. gigas) or ion regulation and Ca transport (C. virginica). To the best of our knowledge, our data provide the first demonstration that hemocytes contain metastable calcium carbonate mineral phases, which are potentially transported to the sites of mineralization, as well as further evidence that only certain fractions of hemocytes are responsible for mineral sequestration and transport. Higher investment in biomineralization function is also apparent from the elevated expression of biomineralization-related genes in other organs (the mantle and the gills) of C. virginica versus C. gigas. These species-specific transcriptomic differences in biomineralizing (mantle and hemocytes) and ionoregulatory (gills) organs are manifested in the oysters' phenotype as superior protective properties of the exoskeleton in C. virginica (shown in this study) or more effective immune defense in C. gigas (Barbosa Solomieu et al., 2015; Elston, 1993; Goedken et al., 2005; Salvi et al., 2014). Further studies are needed to gain insight into the evolutionary pathways and selective pressures that resulted in the species-specific phenotypic specialization on internal (immunity) or external (exoskeleton) protection and assess whether it is based on a true energetic or resource trade-off.

We thank Annett Grüttmüller for performing the routine NanoSIMS operations. The high-resolution TEM studies were carried out at the Nanofabrication and Nanocharacterization Facility, Peterson Institute for Nanoscience and Engineering, University of Pittsburgh, USA.