The gut microbiota of Cystidicola farionis parasitizing the swim bladder of the nosed charr morph Salvelinus malma complex in Lake Kronotskoe (Kamchatka, Russia)

Four individuals of the nosed charr (Fig. 1B) with total weight 370.1 ± 95.5 and total length 345.0 ± 34.2 mm were collected from September 30 to October 01, 2020 in the littoral zone of Lake Kronotskoe (Fig. 1A), which is the biggest lake in the Kamchatka Peninsula (Russia, 54°47′11″ N 160°13′36″ E). The lake area is about 246 km² and has an average depth of 58 m with maximum 136 m. Many rivers and streams flow into the lake, the largest of which are the Listvennichnaya, Unana and Uzon rivers. The lake drains into the Kronotskaya River which flows 39 km southeast into the Pacific Ocean (Busarova et al., 2017). Fish were captured using gill-nets (mesh sizes 35 and 45 mm) and transported alive to the laboratory in plastic containers filled with water from the site of fish capture (duration approximately 30 min). All fish were sacrificed and aseptically dissected as previously described (Kashinskaya et al., 2015). Male and female fish were identified according to gonadal development and the type of morph (N1g or N2) was assigned (Table 1). The level of infestation by C. farionis in the swim bladder of analyzed fish were registered.

Figure 1:

The swim bladder was removed aseptically from the body cavity of fish and placed into a sterile Petri dish. C. farionis individuals (Fig. 1C) were retrieved from the swim bladder using sterile forceps. Before dissection, the nematodes were washed three times in sterile Ringer buffer (pH = 7.4) to remove surface-adherent particles of mucus from the host swim bladder and afterward placed into a sterile Petri dish. In total, 115 individual nematodes were dissected and their guts were removed using a sterile microsurgical forceps, scissors, and scalpel. The “head” end of the nematode was fixed by forceps in order to cut it with scissors. After that, the same manipulations were repeated with the “tail” end of the nematode. Then, the gut was pulled out from the “body” by sterile forceps. All procedures were conducted under stereomicroscope with 8–64× magnification (Carl Zeiss, Stemi^TM DV-4). Finally, “body” remains of dissected nematodes, their guts, and mucosa from the nosed charr swim bladder were taken separately to analyze the associated microbiota. By using the term “body” hereinafter it will be understood this refers to the cuticle, epidermis and longitudinal muscles, etc. which remained after dissection of nematodes and isolation their gut. Taking into account that washing of nematodes in different solution does not mean a complete removal of bacteria from their surfaces and that muscles of nematodes have been considered to be sterile, we can then conclude that the microbiota of the “body” is mostly the microbiota of the external surfaces (cuticle) of the nematode. To check our assumptions, we used micro surgery instruments, as a more relevant method to separate the gut microbiota from microbial communities of the other nematode’s parts (cuticle, epidermis and longitudinal muscles, etc).

The samples of “body” and guts were collected in pools ranging from five to twenty individuals of nematodes (depending on their size) from each studied fish. The number of pooled biological replicates ranged from one to four for each fish although the actual number depended upon the level of infestation. (Table 1). All samples were frozen in liquid nitrogen immediately after dissection and transferred to ‒80 °C storage prior to DNA extraction.

During June–August of 2013, June-November of 2014, and September-November of 2020, 323 individuals of nosed charr were collected in order to estimate the prevalence and intensity of infestation by C. farionis. All fish were dissected and swim bladder was extracted and the number of C. farionis was registered. The prevalence (P), intensity, and mean abundance of parasite infestation were calculated according to the definitions by Bush et al. (1997).

The prevalence (P) of parasite infestation (in %) was calculated as P = I*100/N, where I is the number of infested host and N is the total number of hosts examined. Mean intensity of infestation (I) was assessed as the average of number of individuals of a particular parasite species (K) in a single infested host (n): I = K/n.

For genetic characterization, Cystidicola farionis collected from two nosed charr, were preserved in 96% ethanol and stored at 4°C until DNA extraction. Total DNA was extracted from single ethanol-preserved individuals of parasites using the DNA-sorb B kit following manufacturer’s protocols (kit for DNA extraction, Central Research Institute of Epidemiology, Russia). The 28S nuclear ribosomal large subunit rRNA gene (28S rRNA) was amplified using the following primers (forward 5′-AGCGGAGGAAAAGAAACTAA-3′, reverse 5′-TCGGAAGGAACCAGCTACTA-3′) and PCR conditions as described in Nadler et al. (2006). The fragment of the 18S nuclear ribosomal small subunit rRNA gene (18S rRNA) was amplified using the primers (forward 5′-AGCGGAGGAAAAGAAACTAA-3′, reverse 5′-TCGGAAGGAACCAGCTACTA-3′) and PCR conditions as described in Floyd et al. (2005). Double-stranded DNA was amplified using BioMaster HS-Taq PCR-Color (2x) kit (Biolabmix, Novosibirsk, Russia) according to the manufacturer’s instructions (http://biolabmix.ru/products/klassicheskaja_pcr/biomaster_hs-taq_pcr-color__2_/). The PCR products were purified by adsorption on Agencourt Ampure XP (Beckman Coulter, Indianapolis, IN, USA) columns and subjected to Sanger sequencing using the BigDye Terminator V.3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA, USA) with subsequent unincorporated dye removal by the Sephadex G-50 gel filtration (GE Healthcare, Chicago, IL, USA). The Sanger products were analyzed on an ABI 3130XL Genetic Analyzer (Applied Biosystems). The purification and sequencing of PCR products were performed in the SB RAS Genomics Core Facility (Novosibirsk, Russia). Manual editing, alignment of sequences, distances and phylogenetic analysis were performed with MEGA 7 (Kumar et al., 2016).

Phylogenetic reconstruction was performed using the maximum-likelihood method with the K2 + G model of nucleotide substitutions. Test of phylogeny were calculated using the Bootstrap method with 1000 replications. The nearest-match sequences from GenBank were used for phylogenetic reconstructions (Table 2). Newly obtained sequences were deposited into GenBank (NCBI) under the following accession numbers: MZ151867–MZ151869 for large subunit ribosomal RNA gene, MZ150507–MZ150509 for small subunit ribosomal RNA gene.

Before DNA extraction, parasites, “body” remains and swim bladder were collected into sterile microcentrifuge tubes with lysis buffer (300 µl) for DNA isolation and mechanically homogenized by pestle for 1 min. Following the kit manufacturer protocols, DNA was extracted using the DNA-sorb B kit (NextBio, Russia). The DNA extraction protocol was previously described in Kashinskaya et al. (2020). DNA from a sample containing only sterile deionized water was extracted and included in PCR as a negative control. Before amplification and sequencing of V3–V4 variable region of the 16S ribosomal RNA (rRNA) gene, samples with isolated DNA were mixed in equimolar concentration depending on the type of analyzed sample (Kashinskaya et al., 2015). Sequencing of the V3, V4 hypervariable regions of 16S rRNA genes was carried out on an Illumina MiSeq sequencing platform (500 cycles - 2 × 300 paired-end) by Evrogen (Moscow, Russia) using the primer pair S-D-Bact-0341-b-S-17, 5′-CCTACGGGNGGCWGCAG-3′ and S-D-Bact-0785-a-A-21, 5′-GACTACHVGGGTATCTAATCC-3′ (Klindworth et al., 2013).

The amplification conditions and other methods were applied according to the original manufacturer’s protocol (smetagenomic-libraryµl of reverse primer (1 µM), 5 µl of forward primer (1 µM) and 12.5 µl of 2× KAPA HiFi Hotstart ReadyMix (KAPA Biosystems, Wilmington, MA, USA) in a total volume of 25 µL. The PCR reaction was performed on a 96-well 0.2 ml PCR plate (Life Technologies) using the following program: 95°C for 3 min, followed by 25 cycles of 95°C for 30 s, 55°C for 30 s and 72°C for 30 s and a final extension step at 72°C for 5 min. After producing amplicons, the libraries were cleaned up and mixed in equimolecular portions using SequalPrepTM Normalization Plate Kit (ThermoFisher, Cat # A10510-01) and checked using capillary electrophoresis. Samples were multiplexed using a dual-index approach with the Nextera XT Index kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s instructions.

Forward and reverse read pairs were merged and quality filtered with Mothur 1.31.2 (Schloss et al., 2009). Any reads with ambiguous sites and homopolymers of more than eight bp were removed, as well as sequences shorter that 350 or greater than 500 bp. QIIME 1.9.1 (Caporaso et al., 2010) was used for the further processing of the sequences. De novo (abundance based) chimera detection using USEARCH 6.1 (Edgar, 2010) was applied to identify possible chimeric sequences (‘identify_chimeric_seqs.py’ with an option ‘-m usearch61’ in QIIME). After chimera filtering, the QIIME script ‘pick_open_reference_otus.py’ with default options was used to perform open-reference OTU picking by UCLUST (Edgar, 2010), taxonomy assignment (UCLUST), sequence alignment (PyNAST 1.2.2; Caporaso et al., 2010) and tree-building (FastTree 2.1.3; Price et al., 2010). This algorithm involves several steps of both closed-reference and open-reference OTU picking followed by taxonomy assignment, where the SILVA core reference alignment (release 132; Quast et al., 2013) was used as a reference. Chloroplast, mitochondria and non-bacterial sequences were removed from further analysis. An addition, the singletons were removed from feature table.

The richness (number of OTU’s and Chao1 index) and diversity estimates (Shannon and Simpson index) per sample were calculated using QIIME 1.9.1 (Caporaso et al., 2010). Then, the samples were rarified to the lowest sequencing effort (1667 sequences) using QIIME. Such sequencing effort allowed us to include as many samples as possible in further analyses without reducing the power of the statistical methods. Increasing the sequencing effort did not change the obtained results, but reduced the number of analyzed samples. Nucleotide sequences were deposited in the Sequence Read Archive (SRA NCBI), accession number PRJNA727704.

All data are presented as a mean ± standard error (SE). For estimating the differences between the richness and diversity estimates, a non-parametric Kruskal–Wallis test with Dunn’s multiple comparisons test was applied. A weighted UniFrac (Lozupone and Knight, 2005) dissimilarity matrix was calculated in QIIME and used for downstream analyses. The matrix was used to perform principle coordinates analysis (PCoA) to visualize differences among groups of samples. Permutational multivariate analysis of variance using distance matrices was used as implemented in the ‘adonis’ function of the vegan R package (Oksanen et al., 2018). Pairwise comparisons for all pairs of levels of used factors were performed using ‘adonis.pair’ function of the EcolUtils R package (Salazar, 2018). Analysis of multivariate homogeneity of group dispersions (variances) to test if one or more groups is more variable than the others, was performed using the ‘betadisper’ function of the vegan R-package. In all the aforementioned tests statistical significance was determined by 10 000 permutations.

Parasitological analysis of 323 individuals of nosed charr observed at different years (2013, 2014 and 2020) in Lake Kronotskoe has revealed that the prevalence of C. farionis ranged from 63 to 100%. The intensity of parasite infestation ranged from zero to 550 (Table 3). Prevalence levels observed in different months of 2014 were high and remained stable during all the sampling time.

Four new sequences (two 28S rRNA and two 18S rRNA sequences) were obtained from C. farionis parasitizing the swim bladder of nosed charr. Two sequences (28S rRNA and 18S rRNA) obtained from the C. farionis parasitized in the smallmouth charr were also used in the analysis. The 28S rRNA fragment with length of 866 base pairs were identical to C. farionis from GenBank under accession number MT086834.1. According to the 18S rRNA obtained with sequence length of 756 base pairs a phylogenetic ML-tree was constructed with Ascarophis adioryx as outgroup taxa (Fig. 2).

Figure 2:

Together with the previously published sequences obtained from C. farionis, the original sequences formed a separate clade. The low value of bootstrap support observed for this clade was due to the low variability of the 18S rRNA gene. In general, the C. farionis clade was located on the tree as part of a common weakly resolved clade with Ascarophis adioryx, Comephoronema werestschagini, and the genus Capillospirura.

Eight, nine and ten phyla were registered in the associated microbiota from swim bladder of nosed charr, and the gut and “body” of nematodes, respectively. In all samples the dominant phyla were Proteobacteria for which abundances ranged from 67.9 to 99.7%. The phyla Actinobacteria, Firmicutes and Fusobacteria were also found in a significant amount (13.7, 19.3 and 20.2%, correspondingly) of the associated microbial communities, but their relative abundances were dependent on the individual hosts. Thus, Actinobacteria and Firmicutes of both intestine and “body” were more abundant in the associated microbiota of individuals 2 and 3 in comparison with the other fish. Fusobacteria were more abundant in the associated microbiota of swim bladder than in microbiota of the gut and “body” of nematodes (Fig. 5A). At the lowest taxonomical level examined, the microbiota of gut and “body” of nematodes was mainly represented by Aeromonas, Pseudomonas, Shewanella, and Yersinia. The dominant microbiota of the swim bladder were Acinetobacter, Cetobacterium, Pseudomonas, Shewanella, Yersinia and unclassified Enterobacteriaceae (Fig. 5B). The relative abundances of the main dominant (with abundance more than 3%) of nematode microbiota also varied between individuals of fish. Thus, Cutibacterium, Lawsonella, Staphylococcus, and Tumebacillus were more abundant for both intestine and “body” of nematodes in individuals 2 and 3 in comparison with other fish. The differences in the bacterial community of nematodes parasitizing the swim bladder of nosed charr also depends on the type of tissue analyzed. Thus, Paracoccus were registered only in the microbiota associated with the swim bladder, and bacteria from the genera Lawsonella, Staphylococcus, Tumebacillus and unclassified Alicyclobacillaceae were detected in the microbiota of the intestine and “body” of nematodes, in comparison with microbiota of swim bladder.

Figure 5:

According to the ADONIS test (data not shown) the microbiota associated with the gut of nematodes was not significantly different from microbial communities of their “body” or swim bladder (p > 0.05), nor between different morphs of fish (p = 0.057). But, the principal coordinates analysis (PCoA) also showed a clear grouping of gut and “body” of nematodes in comparison with microbiota associated with the swim bladder of fish, as well as microbiota of nematodes from morph N1g and N2 (Figs. 3 and 4B). The lack of significance of these observations may be due to the strong inter-individual variation in microbiota composition.