C. circumdatus egg masses, larvae, and pupae were collected in November 2018 from plants (Eichhornia crassipes) and riverbank sediments from urban settings surrounding two rivers in Pune, India: (i) the Mula River (18.5551° N, 73.8618° E) and (ii) the Mutha River (18.2901° N, 73.4956° E) (Laviad-Shitrit et al., 2020). Samples were also collected from (iii) laboratory cultures being reared in the Animal House Facility at the Department of Zoology, SPPU, Pune, India.

For the laboratory culture, fourth larval instars were picked from sediments at the Mula riverbank 8 months prior to the sampling for the current study. Larvae collected from the field were reared in the laboratory, and egg masses, larvae, and pupae from subsequent generations were used in this study as laboratory samples (hereafter, laboratory culture). Cultures were maintained in non-toxic plastic tubs (diameter = 35 cm) in an animal house facility at the Department of Zoology (SPPU, India) following ethical guidelines. Culture tubs were layered with sterilized beach sand as an inert substratum material at the base and were provided with fresh water and food every alternate day. Tap water was boiled and filtered to make it potable and microbe-free before feeding the animals. Food mixture comprised of Sphagnum moss and baker’s yeast (5:1). The moss was washed in sterile water, soaked overnight in hot water (70–80°C) and thoroughly rinsed the following day before grinding it to prepare the food diet. Larval transition to the pupal stage was monitored and mature pupae were collected and maintained in culture tubs inside cages. Adults that emerged were allowed to mate and egg ropes laid on the water surface were carefully collected into fresh tubs to raise the progeny.

All the samples that were investigated in the current study were taxonomically identified as C. circumdatus by Laviad-Shitrit et al. (2020). Each sample was treated individually after collection by washing it five times in 1 ml sterile water to remove the microorganisms that were not tightly connected to the samples. All samples were then preserved in storage tubes containing absolute ethanol at –20°C until DNA extraction. Overall, 102 C. circumdatus specimens were collected. The details regarding the number of the specimens and their life stage and sampling location are summarized in Supplementary Table 1.

To remove ethanol residues, samples were centrifuged for 30 min at maximum speed, and the ethanol was removed using a sterile pipette tip. Each individual sample was then incubated in a heat block for 10 min at 100°C, washed in sterile saline water, and crushed with a sterile homogenizer.

DNA was extracted from each sample using a DNA isolation kit (DNeasy Blood and Tissue, Qiagene, Germany) according to the manufacturer’s instructions with minor modifications, as was described previously (Laviad-Shitrit et al., 2017). DNA samples were stored at –20°C.

Genomic DNA was PCR amplified with primers targeting the V4 region of the bacterial 16S rRNA gene. The primers were modified from the primer set employed by the Earth Microbiome Project; CS1_515F (ACACTGACGACATGGTTCTACAGTGCCAGCMGCCGCGG TAA) and CS2_806R (TACGGTAGCAGAGACTTGGTCTGGA CTACHVGGGTWTCTAAT) (Caporaso et al., 2012) (Sigma Aldrich, Israel). The primers contained 5′ common sequence tags (known also as common sequence 1 and 2, CS1 and CS2). Amplicons were created using two-stage “targeted amplicon sequencing (TAS)” as described previously (Naqib et al., 2018).

PCR was performed as described previously (Sela et al., 2020). Sterile DNA-free water was used as a negative control for DNA extraction and PCR amplification to verify that there was no contamination. No contamination was found.

Additional PCR amplification was performed in 10 μl reactions in 96-well plates using MyTaq HS 2X master mix (Bioline, London, United Kingdom). Each sample had a separate primer pair that contained a unique 10-base barcode obtained from the Access Array Barcode Library for Illumina (Fluidigm, South San Francisco, CA, United States; Item# 100–4876). These Access Array primers contained the CS1 and CS2 linkers at the 3′ ends of the oligonucleotides. The conditions for the second PCR and the procedure of the Illumina sequencing were as described in Sela et al. (2020). The second PCR amplification and Illumina MiniSeq sequencing were performed at the Genome Research Core (GRC) at the University of Illinois at Chicago (UIC). Sequence length in the Illumina MiniSeq mid-output flow cell was 2 × 150 paired-end reads.

Overall, 408 fastq files were obtained, corresponding to 102 samples (four files for each sample), with two pair-end sequences each. Data quality of the samples was tested with the fastQC¹ program which is a quality control tool for high throughput sequence data. All the samples were found to be of high quality.

Data pre-processing was performed using the DADA2 pipeline (DADA2 package version 1.14.0, Callahan et al., 2016). A detailed description of the pre-processing steps is given in Laviad-Shitrit et al. (2021). Following data cleaning, both runs were merged by sample, and checked for chimeras. Suspected chimeras were detected and removed using the command “removeBimeraDenovo.” Then a count table for each amplicon sequence variant (ASV) in each sample was produced. Taxonomic assignment of the ASVs was done using DADA2 assignTaxonomy command with the Silva rRNA database (version 138) as reference and minimum bootstrap value of 80%. Then, ASVs of non-bacterial origin (Archaea, chloroplasts, and mitochondria, as well as unclassified taxa at the phylum level), were filtered out. In addition, all ASVs with sequence lengths below 260 bp or above 262 bp were removed. Following this filtering step, sequences were binned into 3,366 bacterial ASVs and then, ASVs with < 50 reads across the dataset were removed from all analyses. 1,959 ASVs were left after this filtration step. The sequencing coverage per sample (read number) after each of the filtration steps is presented in Supplementary Table 2. The distribution of the abundance of each ASV can be found in Supplementary Table 3. Supplementary Table 3 demonstrates the justification for the deletion of ASVs with < 50 reads across the dataset. The final feature table of the ASVs’ taxonomic classification and abundances within each sampling group after the filtration steps (ASVs > 50) is presented in Supplementary Table 4.

Raw reads were recovered as fastq files and are available in the NCBI database under BioProject accession number PRJNA678840.

Rarefaction curves were calculated for each sample using Past 4². The curves demonstrate that the sequencing effort was sufficient to effectively represent the bacterial community composition (Supplementary Figure 1). To determine the effect of sampling location and chironomid developmental stage on microbiota composition, Permutation Analysis of Variance (999 permutations) was conducted using the R software (R Core Team, 2018) vegan package (Oksanen et al., 2007). The model was designed with the developmental stage nested within the sampling location and the bacterial community dissimilarity matrix was calculated using the Bray-Curtis index. Non-metric multidimensional scaling (NMDS) based on Bray-Curtis distances was calculated using PRIMER v.7 (Clarke and Gorley, 2015). NMDS visualizes the differences in the bacterial communities’ composition between the different life stages and sampling locations.

Venn diagrams that represent the distinctive and shared ASVs between each of the developmental stages and locations were created by InteractiVenn software³ (Heberle et al., 2015).

Bacterial richness (Chao1) and diversity (Shannon H’) were calculated using the alpha-diversity analysis tool in the MicrobiomeAnalyst website⁴ (Dhariwal et al., 2017) and visualized with Microsoft excel 2019. These parameters were compared by non-parametric tests (Kruskal-Wallis), followed by post hoc tests using two-sided Mann–Whitney tests adjusted for multiple comparisons with Bonferroni correction (IBM SPSS v.25.0.0.1). The results are presented as mean ± standard error of the mean (SEM).

Linear discriminant analysis Effect Size (LEfSe) was also performed using the microbiome analyst tool (Dhariwal et al., 2017). The analysis was applied on rarefied count data in order to identify the bacterial ASVs that contributed to differences in microbiota composition between the three sampled locations.

The microbiota compositions of the three life stages of C. circumdatus that were sampled from three environments were studied using Illumina sequencing of the 16S rRNA gene.

Similarities between the life stages from the different sampling locations were examined using NMDS (Figure 1). Permutation analysis of variances revealed significant effects of both the chironomid developmental stage and the sampling location on the bacterial community composition (p < 0.001, Table 1). The model attributed approximately 54% of the variation in bacterial community compositions to these two factors. The developmental stage had a more prominent effect [F_{(6, 93)} = 13.76, R² = 0.41, p = 0.001] compared to the sampling location [F_{(2, 93)} = 13.23, R² = 0.13, p = 0.001]. Nevertheless, the bacterial communities of egg masses from the different sampling points clustered together (Figure 1). In addition, the environmental larval samples either clustered together or differed from the bacterial community composition of the laboratory larval culture.

No significant differences were found in microbial richness (Chao1) for the three sampling locations (Mula River, Mutha River, and laboratory culture) (Kruskal-Wallis: H₂ = 0.13, p = 0.94). Similar results were obtained when the bacterial diversity (Shannon H’) was compared for these three sampling locations as well (Kruskal-Wallis: H₂ = 4.45, p = 0.11). However, significant differences were found between the microbial richness (Chao1) of the three life stages (Kruskal-Wallis: H₂ = 46.99, p < 0.001). Richness was the highest in the egg masses stage (2.2, and 1.95 times greater compared to the richness of the larval and pupal stages, respectively). Richness values were significantly different between the egg masses and the larvae, and the egg masses and the pupae (Mann–Whitney: U = 29, p < 0.001; U = 71, p < 0.001, respectively) (Figure 2).

When bacterial diversity (Shannon H’) was compared between the three developmental stages (egg masses, larvae, and pupae), significant differences were observed (Kruskal-Wallis: H₂ = 330.46, p < 0.001) most notably between the egg masses and the larvae (Mann–Whitney: U = 176, p < 0.001) and the egg masses and the pupae (Mann–Whitney: U = 96, p = 0.001). Shannon diversity was the highest in the egg masses (1.8, and 2 times greater than the diversity of the larval, and pupal stages, respectively) (Figure 2).

A total of 1,959 ASVs were detected and of those, 1,026 ASVs were classified into 306 known bacterial genera. The other 933 ASVs (47.62% of the total ASVs) did not have a taxonomic identification at the genus-level. The average prevalence of the genus-level unclassified ASVs was 37.58% out of all the sequences. Aeromonas was the most prevalent genus, detected in 100 samples, with a mean relative abundance of 4.77 ± 1.18%. The most abundant genus was Cetobacterium with 14.20 ± 2.03% relative abundance; however, it was less prevalent than Aeromonas, and it was identified in only 88 samples. Other abundant genera were Dysgonomonas (3.69 ± 0.61%), Vibrio (3.49 ± 1.1%), and Flavobacterium (2.39 ± 0.34%).

Overall, 31 bacterial phyla were detected across the entire dataset. Proteobacteria, Bacteroidetes, and Firmicutes were the most prevalent (detected in all 102 samples) followed by Fusobacteria (98 samples), Epsilonbacteraeota, Actinobacteria (96 samples each), and Cyanobacteria (90 samples). Proteobacteria was also the most abundant phylum with 45.39 ± 2.49% of the ASVs followed by Bacteroidetes (18.81 ± 1.12%), Fusobacteria (14.54 ± 2.01%), and Firmicutes (10.35 ± 1.12%).

Proteobacteria was the most dominant phylum in all the environmental and laboratory egg masses and pupal samples and in the laboratory-reared larvae (Figure 3 and Supplementary Table 5), while Fusobacteria was the most dominant phylum in the environmental larvae (Figure 3). Firmicutes was found in relatively high abundance in the larval samples from Mula River and from the laboratory (17.00 ± 3.12%, and 15.42 ± 1.58% respectively). ASVs belonging to the phylum Epsilonbacteraeota were relatively abundant in the pupal samples from Mutha river (16.40 ± 2.97%) while ASVs belonging to the phylum Cyanobacteria were abundant in the egg masses from Mula River and from the laboratory culture (8.30 ± 3.44%, and 9.46 ± 2.01%) (Figure 3 and Supplementary Table 5).

The genera Flavobacterium and Planktothricoides SR001 were found in high abundance in all the egg masses from the three sampling locations (4.71–10.97% and 1.79–8.04% of the ASVs, respectively) (Table 2). Cetobacterium was found in high abundance in all the larval samples from the three sampling locations (between 13.52 and 47.62% of the ASVs). Dysgonomonas was found in high abundance larvae that were sampled from in Mula River and the laboratory culture (5.84 ± 1.24 and 5.64 ± 0.58%, respectively) (Table 2). Arcobacter and Aeromonas were found in high abundance in all pupal samples from all the three sampling points (5.05–16.18% and 9.33–17.73% of the ASVs, respectively). The genus Vibrio was found in high abundance in the egg masses, larvae, and pupae that were sampled from the Mula River (13.24 ± 4.28, 6.12 ± 4.18, and 6.89 ± 4.64%, respectively) (Table 2).

The distinctive and shared ASVs of the different developmental stages from the different sampling points were analyzed using Venn diagrams (Figure 4). Only a minority of the ASVs were shared between the three life stages in the different sampling locations (Mula River, Mutha River, and laboratory culture: 12.1, 4.1, and 9.7%, respectively, Figure 4A). The highest percentages of shared ASVs were observed in the egg masses from the three sampling points (32.4–67.7% unique ASVs, Figure 4A). Egg masses and pupae from the Mula River shared the highest portion of ASVs (33.4%). In contrast, in the Mutha River, egg masses and larvae shared the highest portion of ASVs (15.2%, Figure 4A). Comparisons of the shared and the unique ASVs of each life stage between the different sampling points showed that life stages from the different sampling points harbored a small portion of shared ASVs (Figure 4B). For each life stage, only a minority of the ASVs were shared across the three sampling locations (egg masses, larvae, and pupae: 17.8, 8.0, and 7.0%, respectively) (Figure 4B).

To identify which ASVs contributed significantly to the variation in the microbiota of the different sampling locations, we used LDA effect size (LEfSe) (Supplementary Figure 2 and Supplementary Table 6). The ten ASVs that most significantly contributed to differences in microbiota composition between the three sampled locations were: Burkholderiaceae (ASV 04), C39 (Rhodocyclaceae, ASV 14), Vibrio (ASV 07), Burkholderiaceae (ASV 37), Arcobacter (ASV 21), Sphaerotilus (ASV 22), Bacteroidia (ASV 12), Flavobacterium (ASV 29), Gottschalkia (ASV 10), and Bacteroidia (ASV 28). In the Mula River, Vibrio, and Gottschalkia were the genera with the highest abundances of all discriminant ASVs. In the Mutha River, unclassified Burkholderiaceae, C39, Arcobacter, Sphaerotilus, and Flavobacterium were the most abundant of all discriminant ASVs. In the laboratory samples, the two unclassified Bacteroidia ASVs 12 and 28, were the most abundant (Supplementary Figure 2 and Supplementary Table 6).