Metatranscriptomics by In Situ RNA Stabilization Directly and Comprehensively Revealed Episymbiotic Microbial Communities of Deep-Sea Squat Lobsters

ABSTRACT Shinkaia crosnieri is an invertebrate that inhabits an area around deep-sea hydrothermal vents in the Okinawa Trough in Japan by harboring episymbiotic microbes as the primary nutrition. To reveal physiology and phylogenetic composition of the active episymbiotic populations, metatranscriptomics is expected to be a powerful approach. However, this has been hindered by substantial perturbation (e.g., RNA degradation) during time-consuming retrieval from the deep sea. Here, we conducted direct metatranscriptomic analysis of S. crosnieri episymbionts by applying in situ RNA stabilization equipment. As expected, we obtained RNA expression profiles that were substantially different from those obtained by conventional metatranscriptomics (i.e., stabilization after retrieval). The episymbiotic community members were dominated by three orders, namely, Thiotrichales, Methylococcales, and Campylobacterales, and the Campylobacterales members were mostly dominated by the Sulfurovum genus. At a finer phylogenetic scale, the episymbiotic communities on different host individuals shared many species, indicating that the episymbionts on each host individual are not descendants of a few founder cells but are horizontally exchanged. Furthermore, our analysis revealed the key metabolisms of the community: two carbon fixation pathways, a formaldehyde assimilation pathway, and utilization of five electron donors (sulfide, thiosulfate, sulfur, methane, and ammonia) and two electron accepters (oxygen and nitrate/nitrite). Importantly, it was suggested that Thiotrichales episymbionts can utilize intercellular sulfur globules even when sulfur compounds are not usable, possibly also in a detached and free-living state. IMPORTANCE Deep-sea hydrothermal vent ecosystems remain mysterious. To depict in detail the enigmatic life of chemosynthetic microbes, which are key primary producers in these ecosystems, metatranscriptomic analysis is expected to be a promising approach. However, this has been hindered by substantial perturbation (e.g., RNA degradation) during time-consuming retrieval from the deep sea. In this study, we conducted direct metatranscriptome analysis of microbial episymbionts of deep-sea squat lobsters (Shinkaia crosnieri) by applying in situ RNA stabilization equipment. Compared to conventional metatranscriptomics (i.e., RNA stabilization after retrieval), our method provided substantially different RNA expression profiles. Moreover, we discovered that S. crosnieri and its episymbiotic microbes constitute complex and resilient ecosystems, where closely related but various episymbionts are stably maintained by horizontal exchange and partly by their sulfur storage ability for survival even when sulfur compounds are not usable, likely also in a detached and free-living state.

ABSTRACT Shinkaia crosnieri is an invertebrate that inhabits an area around deep-sea hydrothermal vents in the Okinawa Trough in Japan by harboring episymbiotic microbes as the primary nutrition. To reveal physiology and phylogenetic composition of the active episymbiotic populations, metatranscriptomics is expected to be a powerful approach. However, this has been hindered by substantial perturbation (e.g., RNA degradation) during time-consuming retrieval from the deep sea. Here, we conducted direct metatranscriptomic analysis of S. crosnieri episymbionts by applying in situ RNA stabilization equipment. As expected, we obtained RNA expression profiles that were substantially different from those obtained by conventional metatranscriptomics (i.e., stabilization after retrieval). The episymbiotic community members were dominated by three orders, namely, Thiotrichales, Methylococcales, and Campylobacterales, and the Campylobacterales members were mostly dominated by the Sulfurovum genus. At a finer phylogenetic scale, the episymbiotic communities on different host individuals shared many species, indicating that the episymbionts on each host individual are not descendants of a few founder cells but are horizontally exchanged. Furthermore, our analysis revealed the key metabolisms of the community: two carbon fixation pathways, a formaldehyde assimilation pathway, and utilization of five electron donors (sulfide, thiosulfate, sulfur, methane, and ammonia) and two electron accepters (oxygen and nitrate/nitrite). Importantly, it was suggested that Thiotrichales episymbionts can utilize intercellular sulfur globules even when sulfur compounds are not usable, possibly also in a detached and free-living state. IMPORTANCE Deep-sea hydrothermal vent ecosystems remain mysterious. To depict in detail the enigmatic life of chemosynthetic microbes, which are key primary producers in these ecosystems, metatranscriptomic analysis is expected to be a promising approach. However, this has been hindered by substantial perturbation (e.g., RNA degradation) during time-consuming retrieval from the deep sea. In this study, we conducted direct metatranscriptome analysis of microbial episymbionts of deep-sea squat lobsters (Shinkaia crosnieri) by applying in situ RNA stabilization equipment. Compared to conventional metatranscriptomics (i.e., RNA stabilization after retrieval), our method provided substantially different RNA expression profiles. Moreover, we discovered that S. crosnieri and its episymbiotic microbes constitute complex and resilient ecosystems, where closely related but various episymbionts are stably maintained by horizontal ex-and metatranscriptomic analyses of the S. crosnieri episymbiotic microbial communities to (i) delineate their taxonomic composition by increasing sensitivity and eliminating biases in PCR cloning, (ii) systematically estimate the effects of the sampling processes from deep-sea hydrothermal vents that are characterized by unusual chemical and physical conditions, and (iii) obtain comprehensive knowledge on their physiology.

RESULTS AND DISCUSSION
Total RNA (rRNA) analysis of S. crosnieri episymbiotic microbial communities. S. crosnieri sampling was conducted at a 982-m-deep hydrothermal vent site in the Iheya North field in the Okinawa Trough in Japan (see Fig. S1 in the supplemental material) in 2015, and episymbiotic microbial community samples subjected to RNA stabilization either in situ or onboard were obtained from the host setae. Total RNA sequencing with Ion PGM System of six in situ and seven onboard RNA-stabilized samples from different host individuals yielded 5,430,057 and 4,847,064 reads, respectively (Table 1). After filtering short and low-quality reads, 619,694 and 695,480 rRNA reads, respectively, were used for taxonomic analyses.
The relative abundances of the rRNA reads confirmed that the orders Thiotrichales, Methylococcales, and Campylobacterales were the major community members. These three major orders showed relative abundances of Ͼ5% in every sample and accounted for 95.3% to 97.6% in total ( Fig. 1). At the genus level, Sulfurovum was particularly dominant and accounted for more than 99% of the Campylobacterales rRNA reads in all samples (see Table S1 in the supplemental material). Notably, Thiotrichales, Methylococcales, and Sulfurovum sequences were consistently found in a previous study that applied a PCR-based method to S. crosnieri episymbiotic microbial communities at the same site in 2007 (20). Therefore, Thiotrichales, Methylococcales, and Campylobacterales (or Sulfurovum) seemed to be stable as major members in S. crosnieri episymbiotic communities at least at this hydrothermal vent site, irrespective of possible PCR biases.
The Thiotrichales, Methylococcales, and Campylobacterales rRNA sequences contained 199, 132, and 337 SILVA identifiers, which were clustered to 82, 44, and 85 operational taxonomic units (OTUs), respectively, at a 95% similarity threshold (Fig. 2). These numbers suggest that each of the three S. crosnieri episymbiont taxa comprised dozens of species instead of a few specifically adapted species. The OTU compositions were highly similar among the different S. crosnieri individuals. The Thiotrichales and Methylococcales episymbionts contained two (OTUs 250 and 180) and one (OTU 17) abundant OTU, respectively. On the other hand, no specific dominant OTUs were observed in the Campylobacterales (Sulfurovum) episymbionts. The saturation of rarefaction curves for all samples suggested that the entire episymbiotic microbial communities were well sampled (Fig. S3a). Collectively, the large numbers of OTUs observed in each sample and the community composition similarities among the samples indicate that the episymbiotic community of each S. crosnieri individual is not com- posed of clonal descendants of a few "founder" cells that settled down after birth or molting of the host but is more complex than previously thought and maintained by horizontal exchange between host individuals around the hydrothermal vent. Episymbiotic microbial community comparisons between the in situ and onboard RNA stabilization methods. While OTU richness and Faith's phylogenetic diversity were not significantly different between the in situ and onboard RNAstabilized samples (P Ͼ 0.05, Mann-Whitney U test; Fig. S3b and c), the in situ RNA-stabilized samples showed significantly larger Shannon diversity (P ϭ 0.018, Mann-Whitney U test; Fig. S3d). Thus, while the OTU sets of the episymbiotic communities were not different between the sampling methods, skewness in their abundances  decreased during the retrieval process from the deep sea. Weighted UniFrac and the principal coordinate analysis (PCoA) revealed that the in situ and onboard RNAstabilized data were significantly different (P ϭ 0.002, permutational multivariate analysis of variance [PERMANOVA], 999 permutations in each test, Fig. S4). It may be noted that the composition of the in situ total RNA 1 sample was exceptionally close to that of the onboard stabilized samples, but this sample was also characterized by the particularly small number of the mapped reads and OTUs, likely because of problems in sequencing (Table 1). We confirmed that exclusion of in situ total RNA 1 sample did not affect the results of the statistical tests of the OTU richness, Faith's phylogenetic diversity, and Shannon diversity.
Comparison of the relative abundances between the in situ and onboard RNA stabilization methods showed relative increases in Thiotrichales (4.0% on average) and decreases in Methylococcales (11.9% on average) and Campylobacterales (6.5% on average) in the onboard stabilization samples, and these shifts were consistently observed regardless of OTUs among different host individuals (Fig. 2). Thus, as expected, the approximately 3-h sampling process without RNA stabilization likely affected the rRNA abundances of the major taxa in a systematic manner. In particular, Thiotrichales OTU 180 (corresponding to the genus Thiothrix) showed a substantial increase under onboard RNA stabilization (Table S2). Notably, Thiothrix was reported to store reduced sulfur compounds as intercellular sulfur granules and utilize them as energy sources by oxidization to sulfate (discussed later) (33)(34)(35).
Metatranscriptomic analysis of S. crosnieri episymbiotic microbial communities. rRNA-depleted RNA sequencing with Ion PGM System of two in situ and two onboard RNA-stabilized samples from different host individuals yielded 23,392,460 reads in total after quality control ( Table 2). Transcriptomic assembly using all reads provided 146,985 contigs (N 50 , 597 bp; longest contig, 22,476 bp), among which 97,348 contigs were predicted to have 109,253 protein-coding sequences (CDSs). We conducted taxonomic and functional annotations of these CDSs, as well as 49,637 contigs on which no CDSs were predicted. Taxonomic annotations were given to 136,313 (85.8%) sequences (i.e., CDSs and contigs), among which 49,014, 27,321, and 26,141 were assigned to Thiotrichales, Methylococcales, and Sulfurovum, respectively. Functional annotations using the KEGG database annotated 78,748 (49.6%) sequences, which are regarded as genes hereinafter. Among those genes, 64,908 (82.4%) were taxonomically annotated to either of the three major taxa.
Transcript abundances were quantified by mapping the rRNA-depleted RNA sequences to the assembled contigs and calculating transcripts per million (TPM) values. Total TPM values of the transcripts that were taxonomically annotated to the three major taxa ranged from 58.8% to 69.6% depending on the samples (Table 2). These values were smaller than the rRNA (total RNA) abundances of the same three taxa, likely because of difficulties of taxonomic annotation of mRNA. The relative transcript abundances among the three major taxa are shown in Fig. 3. Compared to the relative rRNA abundances (Fig. 1), the Sulfurovum episymbionts showed larger mRNA abundances. This might be because the existence of their closely related genomes in the database enabled more sensitive annotation of Sulfurovum mRNA (although Sulfurovum had the  53,394) and/or because the Sulfurovum episymbionts produce more mRNAs than the others. Finally, comparison of the relative abundances between the in situ and onboard RNA stabilization methods showed that the two RNA stabilization methods similarly affected the rRNA and mRNA abundances ( Fig. 1 and 3), i.e., an overall increase in the Thiotrichales profiles but an overall decrease in the Methylococcales and Campylobacterales (Sulfurovum) profiles.
Metatranscriptomic differences between the in situ and onboard RNA stabilization methods. We then compared the transcriptomic profiles of the three major taxa, namely, Thiotrichales, Methylococcales, and Sulfurovum, in depth. Gene-and pathway-level relative transcript abundances were examined based on the functional KEGG Orthology (KO) annotations (Fig. 4, the left and right panels show abundances by genes and pathways, respectively). Those of the genes taxonomically assigned to the three major taxa are also shown in the second to fourth rows in the same figure. As described already, when the onboard RNA stabilization method was adopted, the relative abundances of Thiotrichales genes increased, and those of Methylococcales and Sulfurovum genes decreased. What was most notable here was that these shifts seemed to have been caused by transcriptome-wide changes that coincided with rRNA abundance changes but not with large changes in the abundance of specific genes (e.g., those required for metabolism at hydrothermal vents). It may also be notable that some genes and pathways went against the global trends between the two stabilization methods, either by chance or because of the physiological effects of sampling from the deep sea. However, as far as we investigated, we did not see any genes and pathways that exhibited concerted trends against the transcriptome-wide trends in all three major taxa, even among energy and carbon metabolism-related genes and pathways that were likely systematically affected by deep-sea sampling (Fig. 4). To compare mRNA abundance of each sample by cancelling out those transcriptome-wide effects in each taxon, the expression profiles were also visualized by normalizing per each taxon and sample (Fig. 5). First, the expression profiles of the same taxa from the replicate samples were most similar, showing the robustness of the results. Second, the expression profiles of the same taxa were most similar between the in situ and onboard RNA-stabilized samples when the transcriptome-wide effects were cancelled out. This also confirmed that the sample retrieval from the deep sea did not substantially affect expression of genes of specific pathways. Thus, we concluded that the transcriptomewide shift due to deep-sea sampling was the major factor that affected the RNA expression profiles. Specifically, a reasonable explanation would be that the rRNA and mRNA of Methylococcales and Sulfurovum were globally degraded faster than those of Thiotrichales during deep-sea sampling, possibly because of global attenuation of biological activities by depletion of carbon and energy sources and temperature and hydrostatic pressure changes. We note that previous studies reported gene-specific transcript abundance changes during deep-sea sampling. A study that applied quantitative PCR to analyze the transcript abundance of the Methylococcales pmoA gene, which is a key gene in methane oxidation, reported a 3.7-fold decrease in onboard RNA-stabilized samples of deep-sea S. crosnieri episymbionts compared to in situ RNA-stabilized samples (16). In this study, we observed a 1.82-fold decrease in the TPM values of the pmoA gene, and this could be in line with the 1.51-fold transcriptome-wide abundance decrease of Methylococcales ( Fig. 3 and 4); therefore, the decrease in pmoA transcript abundance might actually be a part of the transcriptome-wide effects of deep-sea sampling. It was also reported that the expression of stress response genes (e.g., recA and hsp90) increases under onboard RNA stabilization (28,30), although these studies dealt with neither S. crosnieri episymbionts nor deep-sea vent microbes. Thus, we analyzed nine genes (recA, dnaK, hsp90, clpB, clpA, clpX, groEL, rpsA, and atpD) for which the expression levels were expected to be increased by environmental stress. In our results, none of the nine genes showed an increase in all three major taxa but collectively followed the transcriptome-wide trends (Fig. 6). A previous study on endosymbionts of Ifremeria nautili, which inhabits deep-sea hydrothermal vents, reported that their 34 h of incubation with few or no electron donors increased the gene expression of the DnaK chaperone system (e.g., dnaK) and the Clp/Hsp100 family of ATP-dependent proteins (e.g., clpA and clpX) (36). However, there are fundamental differences between this study and the previous studies (e.g., 34-h incubation versus 171-min retrieval, endosymbionts versus episymbionts, and taxonomic compositions), and these differences might have led to different transcriptome-wide and gene-specific responses.
Different metatranscriptomic responses among the three major episymbiotic taxa. We found that the rRNA and mRNA of Methylococcales and Sulfurovum were likely rapidly and globally degraded during deep-sea sampling, while those of Thiotrichales FIG 5 mRNA abundance of pathways normalized for each sample and taxon. The heat map shows the total expression abundance of genes in each KEGG pathway, which was scaled by log 10 and z-score normalized for each taxon and sample.
were kept relatively stable. We were then interested in why Thiotrichales episymbionts could keep their rRNA and mRNA more stable than those belonging to the other two major taxa.
We hypothesized that Thiotrichales episymbionts utilized intercellular sulfur globules as an energy source to escape the global attenuation of biological activities. Many sulfur-oxidizing bacteria belonging to Gammaproteobacteria, including Thiotrichales, form intracellular zero-valent sulfur globules as intermediates in sulfur oxidation (35,37). The TusA and DsrAB proteins are known to be involved in the transfer and oxidization, respectively, of sulfur globules (38)(39)(40). In our metatranscriptomic results, we detected the expression of the tusA and dsrAB genes by Thiotrichales episymbionts (TPM values of more than 2,800 and 50, respectively, in all samples). Furthermore, a previous FISH-NanoSIMS analysis revealed that bicarbonate was incorporated into S. crosnieri episymbiotic microbial communities without electron donors and that bicarbonate was enriched in short and thin filamentous cells that resembled Thiotrichales cells (15). Therefore, the Thiotrichales episymbionts were assumed to be able to fix carbon even if electron donors are not supplied using intracellular sulfur globules, and this ability may help them survive even when sulfur compounds are not usable, likely also in a free-living state during horizontal exchange. Notably, the Sulfurovum episymbionts, which are also thioautotrophs but showed a decrease in mRNA abundances, expressed neither the tusA nor dsrAB genes.
Delineation of the metabolic activities of S. crosnieri episymbionts. Finally, we investigated the physiology and ecology of the S. crosnieri episymbiotic community using the in situ RNA-stabilized metatranscriptomic data. Functional annotation using the KEGG database revealed the expression of genes of two carbon fixation pathways, genes of a formaldehyde assimilation pathway, and genes that exploit five electron donors (sulfide [  [NH 4 ϩ ]) and two electron acceptors (oxygen [O 2 ] and nitrate/nitrite [NO 3 Ϫ /NO 2 Ϫ ]). We could not detect the expression of genes related to hydrogen oxidation or nitrogen fixation; these systems were recently observed in symbiotic chemolithotrophic microbes (41,42).
The Methylococcales episymbionts expressed the soxBXY and sqr genes (TPM values of more than 40 in each of the two in situ RNA-stabilized samples), implying that they can also oxidize reduced sulfur compounds. Notably, this possibility of sulfur oxidiza-  tion by the Methylococcales episymbionts of S. crosnieri has not been proposed thus far, whereas metagenomic studies of other microbial ecosystems containing methanotrophs (47,48) and the genome sequences of Methylobacter whittenbury (49) and Methylomarinum vadi (50) indicated that methanotrophs in general can harbor genes for sulfur oxidization.
(ii) Methane metabolism by Methylococcales. Substantial expression of genes for methane metabolism, specifically for complete methane oxidation and formaldehyde (C 1 compound) assimilation, by Methylococcales episymbionts was observed ( Fig. 7 and Data Set S1). These genes included methane monooxygenase (pmoCAB), methanol dehydrogenase (mxaF), those for dehydrogenation of formaldehyde to formate (fae, mtdAB, mch, and fhc), formate dehydrogenase, and hexulose-6-phosphate synthase and hexulose-phosphate isomerase (hps and hpi) in the ribulose monophosphate (RuMP) pathway. Thus, the Methylococcales episymbionts were supposed to be able to oxidize methane to carbon dioxide and assimilate its intermediate, formaldehyde, as a carbon source. On the other hand, only weak expression of six out of nine serine pathway genes (another C 1 compound assimilation pathway) and the ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO) gene was observed. This implied that the Methylococcales episymbionts were mainly composed of type I Gammaproteobacteria methanotrophs, which do not possess RuBisCO and serine pathways, consistent with previous phylogenetic analysis of pmoA genes (16). Nonetheless, it should be noted that type X methanotrophs, which have RuBisCO and serine pathways (51), might also exist in the community. Expression of the hydroxylamine oxidoreductase gene was also observed, suggesting that Methylococcales episymbionts may oxidize ammonia with methane monooxygenase to nitrite via hydroxylamine.
(iii) Respiration and nitrogen metabolism. Expression of genes for both oxygen respiration and aerobic denitrification by all three major taxa was observed. Regarding oxygen respiration, the expression of cytochrome c oxidase cbb3-type subunit I (cco-NOPQ) and cytochrome c oxidase bd-type subunit I (cydAB) genes was detected in all three major taxa (Data Set S1). Regarding denitrification, the expression of the periplasmic nitrate reductase complex (napAB), two cytochrome cd-type nitrite reductases (nirK or nirS), nitric oxide reductase (nor), and nitrous oxide reductase (nosZ) genes was detected in the Thiotrichales and Sulfurovum episymbionts, while expression of only the nirK and norBC genes was detected in Methylococcales episymbionts. Thus, Methylococcales episymbionts seem to generate nitrous oxide (N 2 O), whereas Thiotrichales and Sulfurovum episymbionts seem to generate nitrogen (N 2 ). The substantial expression of genes involved in denitrification suggests the possibility that the three major taxa utilize nitrate and nitrite as an alternative electron acceptor because of their competition for oxygen.
Expression of cytoplasmic nitrite reductase (nirBD or nirA) and nitrogen-assimilating glutamine synthetase (glnA) genes were detected for all three major taxa. Thus, we assume that nitrite in the periplasm would be transported to and assimilated in the cytoplasm.
(iv) Carbon fixation. Expression of the genes of two inorganic carbon fixation pathways among six pathways (43,52) was observed in the Thiotrichales and Sulfurovum episymbionts. The Thiotrichales episymbionts expressed genes encoding the Calvin-Benson-Bassham (CBB) cycle, including the RuBisCO gene. The Sulfurovum episymbionts expressed genes encoding the reductive tricarboxylic acid (rTCA) cycle, including oor, frd, and acl genes. Thus, the CBB and rTCA cycles are likely utilized by the Thiotrichales and Sulfurovum episymbionts, respectively, which is consistent with the physiological characteristics of the cultured strains that belong to Gammaproteobacteria and Epsilonproteobacteria, respectively (43,(53)(54)(55).
Other notable observations would include that Thiotrichales showed relatively high expression of glyoxylate and dicarboxylate metabolism pathway genes and Methylococcales showed relatively high expression of the pentose phosphate pathway genes (Fig. 5).

Conclusion.
In this study, we performed direct metatranscriptome analysis of episymbionts of S. crosnieri by applying the in situ RNA stabilization equipment. To our knowledge, this is the first study that applied and systematically compared the different RNA stabilization approaches to deep-sea symbiotic microbial communities.
Comparison between the in situ and onboard RNA stabilization methods is highly demanded in applying the metatranscriptomic approach to deep-sea microbial communities, which can hardly be directly investigated by laboratory experiments, for three reasons. First, the onboard RNA stabilization method can give biased results that are different from those obtained by in situ RNA stabilization. Second, by doing so, the transcriptome-wide and gene-specific effects of deep-sea sampling can be discriminated. Third, by positively regarding such methodological differences as environmental perturbation, comparative metatranscriptomics may provide new insights into the physiology and ecology of microbial community members, as we did in this study. Last but not least, we note that these three factors would also be considered in metatranscriptomic analyses of other access-limited ecosystems.
The direct metatranscriptome analysis showed that methane utilization and sulfur oxidation are the main energy metabolisms of the S. crosnieri episymbionts. The dominant yet uncultured Thiotrichales, Methylococcales, and Sulfurovum expressed genes that encode almost all the components of the pathways for these chemoautotrophic systems. It may be notable that no expression of hydrogenase genes was detected in this study, while the genome of Sulfurovum sp. strain NBC37-1, which was isolated from an in situ sampler deployed on the same NBC hydrothermal vent site (i.e., not from the S. crosnieri body surface), carries four hydrogenase genes (55,56). Whether the Sulfurovum episymbionts contribute to hydrogen flux will require further studies.
The results of the taxonomic clustering analysis indicated that the composition of the episymbiotic microbial populations is more complex than previously thought, where each of the three taxonomic clades likely contained diverse species. We anticipate that comparative analysis of metagenome-assembled genomes of the S. crosnieri episymbionts powered by long-read single-molecule sequencing technologies will help better understand their true biodiversity, immigration, and evolution within and among different hydrothermal vent habitats in the Okinawa Trough.

MATERIALS AND METHODS
Sample collection and RNA stabilization. S. crosnieri individuals were collected from a single colony at the NBC hydrothermal vent site in the Iheya North field in the Okinawa Trough (27°47.45= N, 126°5 3.80= E, 982-m depth; see Fig. S1 in the supplemental material) during dive 1773 of the remotely operated vehicle (ROV) HyperDolphine on 18 January 2015. The details of the ROV are shown and described in Fig. S2. The seawater temperature of S. crosnieri habitats around the site was 5.7°C on average. All procedures were performed following the manufacturer's instructions and with default settings unless otherwise noted.
In situ RNA stabilization was conducted by following the procedure described by Watsuji et al. (16). Briefly, S. crosnieri individuals were collected in a sample container of the ROV using a suction sampler at the sampling site (Fig. S2). Then, the seawater in the container was replaced with aqueous sulfate salt solution (25 mM sodium citrate, 10 mM EDTA, and 700 g/liter ammonium sulfate; pH 5.2; colored with phenol red for visibility; prepared at 4°C before the ROV dive) by spontaneous diffusion in 9 min (sampling time, 11:46, stabilization time, 11:55 [Japan standard time {JST}]). The aqueous sulfate salt solution, which was proved to inhibit RNase activity (57), was used to provide a large amount of RNA stabilization solution. The S. crosnieri samples were then retrieved by a research vessel. The samples were washed three times using chilled aqueous sulfate salt solution to remove possibly contaminating bacteria that did not firmly attach to the S. crosnieri bodies. The setae covered with episymbionts were dissected from the host individuals and stored in RNAlater (Qiagen) at 4°C overnight and Ϫ80°C subsequently.
For onboard RNA stabilization, S. crosnieri individuals were collected in another sample container using another suction sampler during the same dive. To remove air and to keep the animals alive until retrieval by being kept in cold seawater, the container was prefilled with chilled filtered seawater and kept on Ϫ20°C ice packs before the dive. After collection, the S. crosnieri individuals were retrieved by a research vessel. The samples were washed three times using chilled sterile artificial seawater ( RNA extraction. Total RNA extraction (for 16S rRNA analysis) and subsequent rRNA-depleted RNA preparation (for metatranscriptomic analysis) were conducted. Total RNA extraction was conducted using all 14 seta samples. RNA was extracted using an RNA PowerSoil Total RNA isolation kit (Mo Bio), treated using DNase I (Qiagen) for 10 min at room temperature, and recovered using an RNeasy minikit (Qiagen). The quantity and quality of the RNA were measured using a Qubit RNA HS assay kit with a Qubit 1.0 fluorometer (Invitrogen) and an Agilent RNA 6000 Pico kit with an Agilent 2100 Bioanalyzer (Agilent). We did not observe eukaryotic 18S and 28S rRNA peaks due to host RNA contamination using the Bioanalyzer. rRNA was depleted from four of the total RNA samples that had good quantities and qualities (two each from the in situ and onboard RNA stabilization methods). The samples were treated with the Ribo-Zero rRNA Removal kit for Gram-Negative Bacteria (Epicentre) and purified using RNA Clean & Concentrator (Zymo Research). Using a Bioanalyzer and Qubit, rRNA depletion was confirmed, and the RNA quantity was measured.
The RNA samples were kept at Ϫ80°C until sequencing.
Sequencing and initial quality check. Sequencing of total RNA was conducted using 13 samples that comprised six in situ and seven onboard RNA-stabilized samples (one in situ sample was not used because of a quality problem). The RNA quantities and RNA integrity number (RIN) values ranged from 301.8 to 500.8 ng and from 5.6 to 6.7, respectively. Barcoded cDNA libraries were obtained using the Ion Total RNA-Seq kit v2 and Ion Xpress RNA-Seq Barcode 1-16 kit (Life Technologies). The quantities and library sizes of the cDNA libraries were measured using a Quant-iT dsDNA HS assay kit with a Qubit 2.0 fluorometer and an Agilent High Sensitivity DNA kit with an Agilent 2100 Bioanalyzer, respectively. Each cDNA library was diluted to 50 pM. Then, 5 l of each solution was pooled to prepare two solutions so that each of the pooled solutions contained either the six in situ or seven onboard RNA-stabilized samples. From each of the two pooled solutions, 25 l was loaded onto an Ion 318 v2BC Chip using an Ion PGM Hi-Q Chef kit with an Ion Chef system. Two sequencing runs were conducted using the Ion PGM Hi-Q Sequencing kit for 200-bp libraries with the Ion PGM System (Life Technologies). Sequencer-specific errors were corrected with Karect (58) (-celltypeϭhaploid -matchtypeϭedit). Low-quality sequence regions and short reads were trimmed (-trim_qual_left 20 -trim_qual_right 20) and removed (-min_len 200), respectively, using PRINSEQ v0.20.4 (59).
Sequencing of the four rRNA-depleted RNA samples was conducted as follows. The RNA quantities and RIN values ranged from 316.8 to 436.8 ng and from 5.9 to 6.9, respectively. cDNA libraries were obtained using the Ion Total RNA-Seq kit v2. Quantities and library sizes were measured as described above. Each cDNA library was diluted to 20 pM. From each of the four solutions, 25 l was loaded onto an Ion 318 v2BC Chip using an Ion PGM Hi-Q View Chef kit with an Ion Chef system. Four sequencing runs were conducted using the Ion PGM Hi-Q View Sequencing kit for 200-bp libraries with the Ion PGM System. Sequencer-specific errors were corrected with Karect (-celltypeϭhaploid -matchtypeϭedit). Lowquality sequence regions and short reads were trimmed (-trim_qual_left 26 -trim_qual_right 26) and removed (-min_len 50), respectively, using PRINSEQ. Note that the low-quality filter was set at a stricter level here for subsequent metatranscriptome assembly. rRNA analysis. Reads generated by the two total RNA sequencing runs (i.e., those from 13 samples) were used for 16S rRNA analysis. Taxonomic annotation of each read was performed using BLASTn from BLASTϩ 2.7.1 (60) against SILVA SSURef NR99 123 (61) (http://www.arb-silva.de/) and by retrieving the top hit that satisfied the criteria E value Յ 1EϪ20, alignment length Ն 180 bp, and identity Ն 97%. If the top hit was attributed to multiple SILVA identifiers typically because a conserved region was sequenced, the most frequent identifier from the same stabilization method was chosen based on Occam's razor. The relative abundance of each SILVA identifier was quantified by the assigned read numbers. Then, the SILVA sequences were clustered with a 95% identity threshold using UCLUST v1.2.22q (62,63). Each cluster was designated an operational taxonomic unit (OTU). The abundance of each OTU was calculated as the sum of those of the contained SILVA identifiers. Alpha diversity (OTU richness, Shannon diversity, and Faith's phylogenetic diversity) and beta diversity (weighted UniFrac) were calculated using skbio package v0.5.6 (scikit-bio.org). To calculate Faith's phylogenetic diversity and weighted UniFrac, representative sequences of OTUs selected by UCLUST were aligned using MAFFT v7.470 (64) and phylogenetic trees were reconstructed using FastTree v2.1.11 (65) with the GTRϩCAT model. Principal coordinate analysis (PCoA) was conducted based on weighted UniFrac using the prcomp package in R.
Metatranscriptomic analysis. Reads from rRNA-depleted RNA sequencing were used for metatranscriptomic analysis. Because reference genomes of the episymbionts were not available, all reads generated by the four sequencing runs were merged and fed into the transcript assembly using Trinity v2.3.2 (45) (-min_contig_length 300). For quantification, we used RSEM v 1.2.30 (66) and Bowtie2 v 2.2.9 (67) through the align_and_estimate_abundance.pl script in Trinity (-est_method RSEM -aln_method bowtie2) to map the rRNA-depleted RNA sequencing reads to the assembled transcript contigs and quantify their abundances as transcripts per million (TPM) values.
Protein-coding sequences (CDSs) on each transcript contig were predicted by TransDecoder (http:// transdecoder.github.io/). The taxonomic and functional annotations of the CDSs were conducted using BLASTp (E-value cutoff of 1EϪ5) against the NCBI nonredundant (NCBInr) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (68) (as of 12 March 2018) databases, respectively. Because the transcript contigs could contain partial genes that were missed by the CDS prediction, contigs without predicted CDSs were also queried against the NCBInr and KEGG databases for taxonomic and functional annotations, respectively, using BLASTx (E-value cutoff of 1EϪ5). For the taxonomic annotation using NCBInr, NCBI taxonomic identifiers (TaxIDs) were assigned to the BLAST hits by cross-referencing against the NCBI database (as of 26 March 2019; ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid). Then, the top hit in terms of E value (with bit scores in a tie) among the hits assigned to 1,751 TaxIDs that belonged to the order Thiotrichales or Methylococcales or genus Sulfurovum was retrieved (because Ͼ99% rRNA sequences assigned to the order Campylobacterales were assigned to Sulfurovum at the genus level [ Fig. S1] and Sulfurovum did not belong to Campylobacterales but to unclassified Epsilonproteobacteria in NCBInr [as of 7 March 2018], the genus Sulfurovum was used instead of Campylobacterales for mRNA taxonomic assignment). For the functional annotation using KEGG, the KEGG Orthology (KO) number of the top hit was used if one was assigned.
KEGG orthology, module, and pathway analyses. For each sample, a TPM value of each KO group was calculated as the sum of TPM values of CDSs and contigs that were functionally annotated to that KO. Likewise, the TPM value of each KEGG module and KEGG pathway was calculated as the sum of TPM values of KOs that belonged to each.
Analyses of the transcript abundances of the following genes were performed using the following KOs: pmoA by K10944, recA by K03553, rpsA by K02945, atpD by K02112, hsp90 by K04079, clpB by K03695, clpA by K03694, clpX by K03544, dnaK by K04043, groEL by K04077, and hupL by K06281. Analyses of the transcript abundances of the following functional categories were performed using the following KEGG modules (M numbers) and Ethics approval and consent to participate. The locations for sample collection were not privately owned or protected in any way, and no specific permits were required for the described field studies and sample collection. The field studies did not involve any endangered or protected species.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.