The study was carried out at a single station (Station 18) on the continental shelf off Central Chile (36°27.86′S, 73°10.12′W), 20 km offshore (Figure 1). This sampling station is part of the COPAS Time-Series Program (http://oceanografia.udec.cl/).

Discrete seawater samples were collected monthly at St. 18 at three water-column depths (5, 50 and 85 m) along the O₂ gradient, from July 2010 to May 2011, onboard the RV Kay-Kay II. Temperature, salinity and O₂ profiles were recorded from near the bottom (93 m depth) to the surface using a SeaBird 25 conductivity-temperature depth (CTD) profiler and SeaBird SB43 O₂ sensor. Water samples for DNA, RNA, nutrients, chlorophyll a and flow cytometry were collected using 10-L Niskin bottles.

The water samples for DNA and RNA extraction were immediately pre-filtered onboard, first through a 20-μm mesh and then through a 1.6-μm GF/A filter (47 mm diameter APFA04700, Millipore), using a peristaltic pump (Cole-Parmer 7520-67 Masterflex). For DNA extraction, 10 L of seawater were filtered. The picoplanktonic community was collected onto a 0.2 μm Sterivex filter (SVGP01050, Millipore), which was filled with lysis buffer (50 mM Tris•HCl, 40 mM EDTA, and 0.73 M sucrose). Seawater samples for RNA extraction were filtered during 10 min, after which the filtrate volume was recorded. The picoplanktonic community was collected in a 0.2 μm Durapore filter (25 mm diameter, Millipore). Filters were immediately transferred to 2 mL cryovials containing 300 μL RNAlater (RNA Stabilization Reagent, QIAGEN).

Nutrients were measured in the laboratory, using standard methods (Grasshoff et al., 1999, modified from Strickland and Parson, 1972), ammonium was determined according to the fluorometric method of Holmes et al. (1999), and chlorophyll a with the method described by Arar and Collins (1997). Microbial community abundance was measured by flow cytometry. For this, 1350 μL of seawater was fixed with 150 μL glutaraldehyde (0.1% final concentration) and frozen until analysis. Subsamples of 500–300 μL were processed on a FACSCalibur flow cytometer with an iron-argon laser of 488 nm (Becton Dickinson), as follow: 500 μl of each sample were analyzed for auto-fluorescent cells detection during 5 min. with an average flux of 35 μl/min, using Forward-scattered light and Side-scattered light (FSC-SCC) parameters to estimate the cell size and the FL3 (670/LP nm) parameter for Chlorophyll-a detection and the FL2 (585/42 nm) parameter for phycoerythrin detection. The non-auto fluorescent community was analyzed through SYBR Green staining. 300 μl of sample were stained during 15 min in the darkness and then run during 3 min with a flux of 16 μL/min, FSC-SSC parameters were also used to estimate cellular size and the FL1 parameter (530/30 nm) to visualize the stained cells.

Genomic DNA was extracted using a modified protocol from Galán et al. (2012). Briefly, 2 mg of lyzozyme in 40 μL lysis buffer was added to each filter and cells were lysed for 45 min at 37°C. Following the lysis, 100 μL of proteinase K and 100 μl SDS 20% were added to each filter and the sample was digested at 50°C for 2 h. The lysate was removed, and each filter was incubated with 1 mL of lysis buffer at 55°C for an additional 15 min. The total pooled lysates were then extracted, once with phenol-chloroform-isoamyl alcohol (25:24:1; pH 8) and once with chloroform-isoamyl alcohol (24:1), and samples were concentrated by isopropanol precipitation.

The RNA extraction was carried out according to a modification of the mirVana miRNA Isolation kit (Ambion®, Life technologies). Samples were thawed on ice and RNAlater reagent was removed and discarded. The cell lysis was performed by adding 100 μL zirconium beads (OPS diagnostics, LLC) and 750 μL Lysis/Binding buffer. After processing the samples in a bead beater, total RNA was extracted according to the manufacturer's protocol. This briefly consists in phase separation with acid-phenol-chloroform, the RNA is then washed within a column, from which is separated from the sample and eluted with elution solution (Ambion®, Life technologies).

RNA was purified from any DNA contamination using the TURBO DNA-free kit (Ambion®, Life technologies). It consists of incubation with DNAse TURBO enzyme at 37°C for 30 min and further addition of DNA inactivation reagent. Possible remaining DNA was monitored by PCR amplification of RNA aliquots. cDNA synthesis was carried out according to the Sensiscript Reverse Transcription kit (Qiagen). Initial RNA concentration was standardized to 20 ng/μL, by dilution with RNAses-free water, to a final volume of 12 μL. After the dilutions were set, the manufacturers protocol was followed, with the preparation of two master mixes, one containing the Sensiscript enzyme, and the other without it (control). The reaction consisted of 8 μL of master mix and 12 μL of RNA template. The reverse transcription was performed with incubation at 37°C for 1 h.

DNA, RNA, and cDNA concentrations were quantified in a NanoDrop ND 1000 spectrophotometer.

A pair of degenerate primers (cbbM-R and cbbM-F; Table 1) was designed to assess the diversity of cbbM genes of putative GSO off central Chile present during two contrasting periods. Primer design was carried out by the alignment of cbbM sequences of ESP GSO provided by previous metagenomic data: Contig14992, (GenBank accession number KJ783409, Murillo et al., 2014), Candidatus Ruthia magnifica (CP000488), Candidatus Vesicomyosocius okutanii (AP009247) and uncultured SUP05 bacterium from Saanich Inlet (ACX30513). Primers were designed with the PRIMER 3 software (v. 0.4.0) and synthesized at IDT (Integrated DNA Technologies, USA).

PCR mixtures were prepared in a total volume of 25 μL containing 14.25 μL of nuclease free water, 5 μL of 5x green GoTaq Flexi buffer, 2 μL of 25 mM MgCL₂, 0.5 μL of 10 mM dNTP, 1 μL of 10 μM forward primer, 1 μL of 10 μM reverse primer, 0.25 μL of 5 u/μL GoTaq® DNA Polymerase (Promega) and 1 μL of template. Optimal thermal conditions for PCR amplification were as follows: initial denaturation at 95°C for 3 min, followed by 11 cycles with an annealing temperature starting at 65°C and decreasing by 0.5°C per cycle, reaching 60°C. Finally, 19 cycles of annealing at 60°C was performed. Denaturation and elongation remained at 95 and 72°C, for 30 s and 1 min, respectively, in every cycle and a final elongation step was performed at 72°C for 5 min. PCR products were visualized on a 2% agarose gel, from which DNA was then purified (Nucleic Acid and Protein Purification kit, Macherey-Nagel) and cloned, according to the pGEM®-T easy Vector System kit (Promega). Two samples were selected, for contrasting environmental periods (March through July, representing the transition from austral summer to winter conditions). Cloning was confirmed by PCR amplification with pucM13F and pucM13R primers. Clones were sent for sequencing at Macrogen Inc. (Geumcheon-qu, Seoul, Rep. of South Korea).

In order to evaluate the cbbM diversity at the study site, a phylogenetic tree was built using 75 clone sequences obtained using degenerated primers and grouped by similarity at 98% generating 8 OTU groups (GenBank: KM275794-KM275837), these cbbM sequences were then, compared with the cbbM GSO from Contig1499281 and 81 sequences comprising the four forms (I, II, III, and IV-like) of RubisCO, plus 45 sequences from the GOS expedition described by Tabita et al. (2007b) and Tabita et al. (2008), together with 42 selected sequences of the form II that cover the diversity of chemolithoautotrophic representatives (Table S2). The amino acid sequences were aligned using MUSCLE 3.6. (Edgar, 2004) and the phylogenetic inference was ran in the software tool Bosque (Ramírez-Flandes and Ulloa, 2008) using maximum likelihood implemented in PhyML (Guindon and Gascuel, 2003) and the WAG substitution model (Whelan and Goldman, 2001).

A pair of specific primers (GSOcbbM-R and GSOcbbM-F; Table 1) was designed for the amplification and quantification of the RuBisCO form II (cbbM) gene of the GSO group off central Chile. The design was carried out by the alignment of the cbbM GSO sequence, encoded in Contig14992 and the sequences of the clones obtained with the degenerated primers. Primers were designed in PRIMER 3 software (v. 0.4.0) and synthesized at IDT.

PCR mixtures were prepared in a total volume of 25 μL, as described above. PCR was performed according to the following thermal cycling: initial denaturation at 94°C, for 5 min, followed by 29 cycles of denaturation at 94°C for 1 min, annealing at 61.6°C for 1 min, and final elongation at 72°C, for 1 min and 30 s. PCR products were always visualized on an agarose gel (1%) stained with GelRed™. A clone library (promega) was performed to test primer specificity; randomly selected clones of the PCR amplicon (n = 100) were sequenced to verify that the primers displayed specificity for cbbM. Seventy five clones contained the targeted sequence, the remaining 15 sequences were divide in 12 low-quality sequences pulled out of the analysis and 3 sequences with low identity to cbbM, low bit-score and high E-value, that were also excluded from the analysis.

Quantitative PCR (qPCR) was performed to assess cbbM abundance of the putative GSO, using DNA as a target. Also, cbbM transcripts abundance was evaluated using cDNA as a target.

Standard curve preparation was carried out through standard dilution series, using a target-specific amplicon, obtained from clone libraries of the sample MLC2285m (1st March, 85 m). For this, the concentration of the PCR product was determined in a NanoDrop ND 1000 spectrophotometer and transformed to numbers of copies per μL. A dilution of the PCR product was preformed to 4 × 10⁷ copies/μL, corresponding to the stock solution, from which a serial dilution is carried out (dilution factor 1:10), with a total of seven dilutions (4 × 10⁷, 4 × 10⁶, 4 × 10⁵, 4 × 10⁴, 4 × 10³, 400, and 40 copies/μL). Negative control corresponds to nuclease free water.

Quantitative PCR was conducted in a StepOnePlus thermocycler (Applied Biosystems) and amplification was detected through the SYBR Green I method. Quantitative PCR mixtures were prepared in a total volume of 20 μL containing 1X of buffer with polymerase enzyme, 0.8 μM of GSOcbbM reverse primer, 0.8 μM of GSOcbbM forward primer, 5.8 μL of nuclease free water and 1 μL of 15 ng/μL template. Every sample were analyzed in triplicate. qPCR was performed according to the following thermal cycling: initial denaturation at 95°C, for 10 min, followed by 30 cycles of denaturation at 94°C for 1 min, annealing at 61.6°C for 1 min, and final elongation at 72°C, for 1 min and 30 s. Standard curve for DNA qPCR showed a slope of −3,5401x + 31,239 along with a R² = 0.997. All values were estimated above a threshold of ΔRn = 0.021. Standard curve for cDNA qPCR showed a slope of −3,6059x + 35,581 with a R² = 0.998. All values were estimated above a threshold of ΔRn = 0.400 (Figure S3).

Final data were exported from the thermocycler into the StepOne Software 2.2, as number of copies per μL. Then, the data were normalized by the volume of filtered water, and presented finally as number of copies per mL.

To establish the relationship between the factors of seasonality and depth with cbbM abundance and transcription, two-ways analyses of variance was performed (ANOVA) with significance defined as p < 0.05. Spearman correlations were calculated to determine the relationship between biological variables (cbbM gene and cbbM transcript abundance) and environmental variables (temperature, O₂, NO3−, NO2−, NH4+, total bacterial abundance and Chlorophyll a), with significance defined as p < 0.05. Analyses were performed in STATISTICA v.10 Software.

Metatranscriptomic datasets for Station 18 were retrieved from iMicrobe (formerly CAMERA) under accession number CAM_P_0000692. The sequences where aligned against the Kegg protein database using Diamond (Buchfink et al., 2015). The Kegg protein database was constructed from a subset of protein sequences with assigned Kegg Orthology (KO). To count enzymes given the EC number, all the corresponding KOs were considered, using the same Kegg database, which contains tables with assignation of Kegg protein sequences to KO and then EC, when it is the case. The relative expression for each gene analyzed is expressed as; 100^*(Ni/Nt) = X, with {X|0.0 < X < 100}. Where Ni is the number of assigned sequences to a given gene and Nt is total number of expressed sequences for the given metatranscriptome.

The study area is characterized by seasonal coastal upwelling and the generation of seasonal O₂ minimum and anoxia. This can be observed in the vertical distribution of hydrographic conditions along the study period (O₂, NO3−, NO2−, bacterial abundance [BA] and NH4+; Figure 2). Vertical distribution of O₂ concentration (Figure 2A) displays a marked seasonality, and during austral winter (June 21st to September 21st), O₂ dropped to <50 μM at 50 m depth, reaching a minimum value of 18 μM at 85 m. Meanwhile, during austral summer (December 21st to March 21st), the oxycline is shallower, ~15–20 m depth, and O₂ concentration dropped sharply below 30 m, reaching undetectable levels at deeper depths. During these months, the microbial abundance showed their maximum values along the water column, especially in the subsurface layers and bottom waters, as measured by flow cytometry (see Section Materials and Methods, Figure 2B), with a maximum of 4.42 × 10⁶ cells^*mL⁻¹ at 5 m, 1.54 × 10⁶ cells^*mL⁻¹ at 50 m, and 1.84 × 10⁶ cells^*mL⁻¹ at 85 m.

In this hydrographic context, the nitrogen species nitrate (NO3−), nitrite (NO2−) and ammonium (NH4+) also show a seasonal variability with a marked decline in nitrate concentration during the upwelling period, which is more pronounced at 85 m, with a concomitant accumulation of nitrite at bottom depths (Figures 2C,D). The reduction of nitrate to nitrite (the first step in denitrification) is coherent with the utilization of nitrate as alternative electron acceptor, with the accumulation of nitrite being a signature of functional anoxia (Thamdrup et al., 2012). The accumulation of nitrite co-occurred with an accumulation of ammonium at 85 m, but also at intermediate depths (Figure 2E), a relationship that was supported by a high correlation value (r = 0.735 p < 0.05, Table S3), consistent with an active nitrogen cycling in these low-oxygen ocean areas (Lam and Kuypers, 2011).

With the set of primers designed to analyze the cbbM diversity, we were able to recover 36 and 39 bona fide clone sequences from 50 to 85 m, respectively. These 75 sequences were grouped by similarity at 98%, and 8 OTU groups were generated and used for further analysis. A phylogenetic inference was constructed comparing our sequences along with 81 sequences for the four forms (I, II, III, and IV-like) of RuBisCO and 45 sequences from the GOS expedition described by Tabita et al. (2007b) and Tabita et al. (2008), together with 42 selected sequences of the form II that cover the diversity of chemolithoautotrophic representatives. Consistently, the 8 OTUs were associated with the RuBisCO form II sequences and 7 of them clustered together, along with the sequence cbbM GSO from Contig14992 (ESP_GSO), described in Murillo et al. (2014) (Figure 3B; bootstrap of 100%). This cluster was named ESSP GSO-St18 (Eastern Subtropical South Pacific GSO—Station. 18) and was significantly separated (bootstrap = 99%) from its closest relatives, the sulfur-oxidizing γ-proteobacteria clade known as SUP05/ARCTIC96BD-19 (Walsh et al., 2009; Canfield et al., 2010; Anantharaman et al., 2013; Glaubitz et al., 2013; Marshall and Morris, 2013), formed by the sequences of SUP05 from Saanich Inlet, symbiont species Candidatus Vesicomyosocius okutanii and Candidatus Ruthia magnifica, Candidatus Thioglobus autotrophicus, four uncultured bacteria and two members of ARCTIC96BD-19 subclade. The remaining OTU was associated to the ARCTIC96BD-19 sequences.

Analysis by qPCR amplification of cbbM gene of the ESSP-ST18 cluster (Figure 4), showed clear evidence of higher abundance during the austral spring-summer, especially at 85 m, with a maximum of 199 ± 19 × 10³ copies mL⁻¹ in December. Conversely, the minimum values of cbbM abundance, at 85 m, were present in July (116.7 ± 39.7 copies mL⁻¹). On the other hand, at 5 and 50 m, there was less variation during the annual cycle. At 50 m, despite the low fluctuation, there was a clear maximum in summer, with 67 ± 28 × 10³ copies per mL. Although the lower variation and the minimum values (5.6 ± 2.5 copies mL⁻¹ in July) were observed at 5 m, the values at 5 m were at least two orders of magnitude less than at 50 and 85 m of depth in all samples.

An ANOVA test indicated that the differences between both depths and months were significant (Table S1, F = 8.432; p = 0.002; F = 2.380; p = 0.047, respectively). Spearman correlation analysis showed a significant negative relationship between cbbM abundance and O₂ concentration (Table S3, r = −0.804; p < 0.05). Thus, cbbM abundance fluctuates during the annual cycle, increasing with the decrease in O₂ concentrations.

The overall microbial community abundance was also maximal during the upwelling season (austral spring-summer; Figure 5). With the intention of exploring the importance of the ESSP-ST18 cluster possessing the cbbM gene in the total microbial community, we normalized cbbM abundance with the microbial community abundance measured by flow cytometry (Figure 6), assuming one copy of cbbM gene per cell of ESSP-ST18 members (Tourova and Spiridonova, 2009). Then we used this normalization as a proxy to estimate the relative contribution of this group to the microbial community, finding the highest percentages (10 to 16%) between December and March, at 85 meters of depth. By contrast, the lower percentage of contribution was in July, at 5 m.

To determine the transcription levels of cbbM of the ESSP-ST18 cluster, qPCR was performed with cDNA as a target. Results (Figure 7) showed a considerable peak during the austral summer, and maximum values at 85 m (20.1 ± 3.1 copies mL⁻¹ in December). However, cbbM transcript values were two to three orders of magnitude lower than cbbM gene abundance. Although the patterns were similar for both variables (Figures 4, 7), it was clear that they were not highly correlated (Table S3, r = 0.440; p < 0.05). Thus, the RuBisCO enzyme did not appear to be transcribed by every member of the ESSP-ST18 group. A maximum value at 5 m during March survey (12.9 ± 1.39 copies per mL) is highlighted and cannot be explained by environmental factors or biological variables, being probably an effect of contamination or manipulation error. As such, it was decided to keep it outside of the statistical analyses.

An ANOVA test indicated that the differences between months were significant (Table S2, F = 6.133; p = 0.001), while the differences among depths were not (Table S2, F = 1.921; p = 0.183). However, cbbM transcript abundance showed significant differences between depths during austral spring-summer months (November, December and January; F = 9.609; p = 0.029). Spearman correlation analysis did not show an important relation between cbbM transcript abundance and O₂ concentration (Table S3, r = −0.231; p > 0.05) or any environmental variables, suggesting a non-linear relationship between cbbM transcript and the variables analyzed in this study. However, there was a clear seasonal pattern in its distribution (Figure 7).

Within the annual cycle explored, we examined six metatranscriptomes for the same depths analyzed above, and for two contrasting periods: three corresponded to downwelling (August 2010; 5, 50 and 85 m) and three to upwelling (March 2011; 5, 50, 85 m). These data allowed us to assess the community expression of key genes, not only for carbon fixation, but also for the nitrogen and sulfur cycles, and to explore the coupling of these cycles by chemolithoautotrophs in this seasonal OMZ. We analyze the percentage of relative expression of the key functions described for the six pathways of inorganic carbon fixation known to date (see Hügler and Sievert, 2011, for a review) in each metatranscriptome and compared them among depths and periods (Figure 8). From all the functions evaluated, only three appeared with a relevant relative expression (>0.1), corresponding to: the enzyme Ribulose 1,5-bisphosphate carboxylase from the CBB cycle and the enzymes Pyruvate synthase and Pyruvate:water dikinase, which participate in several cellular pathways, but together with the enzyme PEP carboxylase, have been related to the dicarboxylate/4-hydroxybutyrate (DC/4-HB) cycle when they act in a concerted way to convert acetyl-CoA to oxaloacetate. Since sequence similarity based inference was through KEGG Orthology (KO), the sequences were associated to the best hit against a KO identifiers (K numbers) and cannot distinguish, in the case of RuBisCO, among the different forms of the large subunit (I, II, or III). Therefore, the results presented in Figure 8 for RuBisCO relative expression are a summation of these forms. Considering the above information, the high levels of expression for RuBisCO at 5 m for both upwelling and downwelling periods are related to the form I of the enzyme and their coding gene cbbL (or also rbcL), associated to a great diversity of photoautotrophs, while at 85 m during upwelling the 82% of the sequences are associated to the Gamma sulfur-oxidizers (GSO) group (Figure 9). The GSO group considered in this work for the assignation of sequences is formed by: SUP05, the symbiont species Candidatus Vesicomyosocius okutanii and Candidatus Ruthia magnifica, Candidatus Thioglobus autotrophicus, Candidatus Thiglobus singularis, members of ARCTIC96BD-19 and thiotrophic symbionts of diverse hydrothermal vent and deep fauna (Table S4). The above was confirmed when extracting the expressed sequences associated to the GSO group for cbbM, rpoC (RNA polymerase subunit beta) and groEL (molecular chaperone), the latter two used as internal controls of relative expression, which showed that the expression of cbbM genes related to GSO group was only relevant at the anoxic depth (Figure 8 inset).

In the case of sulfur and nitrogen cycling coupled to carbon fixation, we evaluate the in situ expression of key genes for chemolithoautotrophic metabolism proposed for the GSO group in the study area by Murillo et al. (2014). Figure 10 shows the percentage of relative expression of the different steps for the oxidation of reduced-sulfur compounds, with nitrate as electron acceptor for the generation of reducing power and ATP. All the pathways related to sulfur cycling showed the maximum expression at 85 m in the upwelling season (March 2011, Figures 10A,B), being the pathways for utilization of elemental sulfur (DsrAB) and sulfite (AprA) the most highly expressed and one order of magnitude higher than the thiosulfate oxidation (sox) genes (Figure 10A). As for RuBisCO, for each pathway we selected a key gene and analyzed the taxonomic composition of the sequences associated along with the percentage of relative expression for the three depths under study (Tables 2, 3); for the sulfur pathways the transcripts were dominated at 85 and 50 m by sequences associated to the GSO group of bacteria (Table 2).

Regarding nitrogen cycling, high transcription levels of the nitrate reductase gene cluster (narGHIJ, Figure 10C) were observed, coherent with the utilization of nitrate as alternative electron acceptor in these low-oxygen, even anoxic, areas of the ocean. Although the expression is still dominated by the GSO group at 85 m, with 62.65% of the sequences, at 50 m takes second place with 23.17% of the expressed sequences, being overtaken by Nitrospina sp. with 46% of the expressed sequences (Figure S1). This result shows that there is a great diversity of organisms carrying out this central pathway at the OMZs. Finally, the expression of functions for the utilization of nitrite shows that nirK was the only relevant gene being observed to be expressed (Figure 10D). However, it did not follow the expression pattern of the other genes, having low expression levels at 85 m during upwelling and the highest relative expression at 50 and 80 m during downwelling, even two orders of magnitude higher than the other pathways analyzed (Table 3). From the few sequences expressed for nirK at 85 m during upwelling, 7.67% were associated to GSO organisms, 9.21% to Thaumarchaea, 10.74% to Nitrospina sp. SCGC_AAA799_C22, 10.74% to an uncultured bacterium, 9.46% to the Alphaproteobacterium Phenylobacterium zucineu, and 3.32% to Nitrosospira sp. NpAV (Figure S2). The rest of the sequences were divided among a wide variety of bacterial genera (data not shown). For the case of the 50 m metatrascriptome during upwelling, and the 85 and 50 m metatranscriptomes during downwelling, archaea dominated the community expression of nirK gene, with a 79.19, 99.79, and 99.82%, respectively. The expression of the nirBD operon was quite low, with only the nirB gene being detected (Figure 10D), but the sequences expressed were dominated by the GSO group at 85 and 50 m with a 79.69 and 84%, respectively. Meanwhile, the genes encoding the enzymes for nitrous oxide (N₂O) production norB and molecular nitrogen (N₂) nosZ showed also a very low expression level, but for both genes the sequences expressed showed only a 4.0 and 4.69% associated to the GSO group. Figure 10E shows a summary of proposed reactions for the coupling of carbon, sulfur and nitrogen cycles at the OMZ by chemolithoautotrophs, with the percentage of relative expression and taxonomic association to the GSO group for each key gene analyzed in this work (figure modified from Murillo et al., 2014).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.