Environmental Factors and Pollution Stresses Select Bacterial Populations in Association With Protists

Digestion-resistant bacteria (DRB) refer to the ecological bacterial group that can be ingested, but not digested by protistan grazers, thus forming a specific type of bacteria-protist association. To test the hypothesis that the environment affects the assembly of DRB in protists, a mixotrophic ciliate, Paramecium bursaria, and a heterotrophic ciliate, Euplotes vannus, were reared at different temperatures, light conditions, and concentration gradients of antibiotic oxytetracycline and heavy metals. Community profiling indicated that the composition of DRB in both species varied significantly across the manipulated conditions, except for in P. bursaria under light/dark treatments. Clone library analysis of bacterial 16S rRNA genes showed that DRB were diverse. Pseudomonas became more abundant during the warmer treatment of P. bursaria, whereas the dominance of Pseudoalteromonas weakened and Vibrio became more abundant in E. vannus at a higher temperature. During the treatment of diel light:dark cycles, Aestuariibacter and Alteromonas were selected for in E. vannus but not Pseudoalteromonas, which was highly represented in the all-light and all-dark treatments. In contrast, P. bursaria consistently hosted Nevskia, Curvibacter, and Asticcacaulis under all light conditions. There were many bacterial species co-resistant to oxytetracycline and to protistan digestion, in which Sphingomonas, Alteromonas, Aestuariibacter, Puniceicoccaceae (Verrucomicrobia), Pseudomonas, and Sulfitobacter were frequently abundant. Flectobacillus and Aestuariibacter were major lead-resistant bacteria associated with the studied protists. Acinetobacter and Hydrogenophaga were abundant in the P. bursaria treated with a high dose of mercury. Aestuariibacter was found as a dominant group of DRB in E. vannus across all cadmium treatments. In summary, this study demonstrates for the first time that environmental stress selects for bacterial populations associated with protists and that there are diverse bacterial species that not only are resistant to pollution stresses but can also survive protistan predation. This work highlights that bacteria-protists associations need to be taken into account in understanding ecological and environmental issues, such as resilience of bacterial community and function, microbial co-occurrence, and quantity and distribution of antibiotic resistant bacteria and genes.


INTRODUCTION
Bacteria and protists are major players in the microbial food web of aquatic environments, in which heterotrophic bacteria incorporate dissolved organic matter and nutrients for growth and transfer carbon and energy to higher trophic levels as a consequence of predation by small protists (Azam et al., 1983). In the last two decades, much efforts have been made to understand how and why protists selectively graze on bacteria or how bacteria become resistant to protistan predation. It has been demonstrated that changes in bacterial cell size, morphology, movement, cellcell configuration in space, and chemical resistance contribute to their anti-predation (Matz and Kjelleberg, 2005;Pernthaler, 2005;Jousset, 2012). Furthermore, recent studies have provided evidence for another strategy of anti-predation-some bacterial populations in aquatic environments can be ingested, but not digested by protistan grazers, thus potentially forming microbial associations with the protistan hosts (e.g., First et al., 2012;Šimek et al., 2013;Gong et al., 2014). When protistan isolates and assemblages are exposed in natural bacterial communities, the detected digestion-resistant bacteria (DRB) can be prevalent and diverse, showing a common trait in taxonomic composition where Gammaproteobacteria and Alphaproteobacteria usually dominate in marine samples (Martinez-Garcia et al., 2012;Pucciarelli et al., 2015;Farnelid et al., 2016;Gong et al., 2016). From an ecological point of view, the existence of DRB could prevent organic matter and energy from being transferred to protistan grazers and higher trophic levels via the microbial loop, adding another layer of complexity in the structure and functions of microbial webs.
Despite the ecological importance of DRB, associations between DRB and protists remain poorly understood. One issue is how environment and host affect the diversity and assemblage composition of the DRB associated with protistan populations. There are indications that pre-culture conditions would affect the ingestion and digestion of heterotrophic nanoflagellates (Boenigk et al., 2001). In our previous studies of DRB in ciliated protists cultivated in the laboratory or isolated directly from natural environments, we also observed that the DRB assemblage could substantially differ among host species and even for populations of the same species isolated at separate time points at an identical location, suggesting that the DRB diversity and assemblage composition can be both host-and environmentdependent (Gong et al., 2014(Gong et al., , 2016. These relationships seem straightforward, since numerous studies have demonstrated that changes in physical and chemical factors drive bacterial abundance and community composition in bulk environmental samples via bottom-up controls or due to effects of pollutants (e.g., Storesund et al., 2015;Ibekwe et al., 2016).
It is well known that environmental factors play important roles in shaping the diversity and dynamic of the bacteria in a natural aquatic environment (Alonso-Sáez et al., 2006;Adams et al., 2010;Darriba et al., 2012;Hu et al., 2014). When the diversity of source bacteria shifts via environmental filtering, it is likely that different sets of bacterial populations are ingested by the protists and retained inside the cells for at least a short-term period, resulting in different assemblage structures of DRB accordingly. However, when the ingested bacteria escape from lysosomes, they enter into the cytoplasm of host cells, a relatively stable physical and chemical environment (due to cellular homeostasis), which may select for a specific group of DRB. However, it remains to be investigated if and how DRB assemblage associated with bacterivorous protists varies with environmental changes and under chemical stresses.
In this study, we performed microcosm experiments using two bacterivorous protist species, Euplotes vannus, a marine heterotrophic ciliate, and Paramecium bursaria, a freshwater ciliate with endosymbiotic green algae (Chlorella sp.). We aimed to test the following hypotheses: (1) warming, light condition, antibiotic, and heavy metal treatments would individually lead to significant changes in DRB assemblages of both ciliate species; (2) the changes in DRB assemblages in both ciliates were greater than the free-living bacterioplankton in the media; and (3) there were diverse bacteria co-resistant to antibiotics/metals and protistan digestion. We believe that the exploration of the diversity and variation pattern of DRB will improve our understanding of bacteria-protist interactions and the dynamics of microbial food webs.

Source Organisms and Cultivation
The freshwater ciliate species P. bursaria, collected from a pond on the campus of Yantai University (Yantai, Shandong, China) in 2014, and the marine species E. vannus, isolated from the estuarine water of the Guangdang River (Yantai, China) in 2013, were investigated. The water used for cultivation was collected on site and treated by passing through a GF/C glass fiber filter (1822-047 grade, pore size: 1.2 µm, Whatman, United Kingdom), such that most pico-, nano-, and larger-sized organisms were removed, and the raw bacterioplankton was largely retained as food for the testing protists. The ciliate cells were observed under a stereoscope (XLT-400, Guilin Guiguang, Guangxi, China), individually selected with a micropipette, and washed for at least 3 times in sterile water to remove other protists and attaching bacteria. A few (3 to 5) cleaned cells were then inoculated into filtered water, allowing the stock cultures to develop. The stock cultures (each in 15 mL volume) were maintained in glass Petri dishes at 16 • C for 2 weeks, with several rice grains to enrich bacterial prey (Supplementary Figure S1).

Experimental Settings
Different temperatures, illuminations, and stresses of antibiotics and heavy metals were set up to investigate the effects of environmental factors and pollution stresses on the diversity of DRB. For each treatment, about 30 ciliate individuals were picked up from the stock cultures and transferred into a Petri dish containing 20 ml aliquote of the filtered water, in which we assumed there was a bacterial community initially identical in quantity and community composition but subsequently being selectively grazed under various treatment conditions (Supplementary Figure S1). Effects of each factor were investigated independently, and each treatment was performed in triplicate. To investigate the effect of temperature, the Petri dishes with stock ciliate cultures were incubated in three chemostats, which were set at a temperature of 16 • C, 21 • C, and 25 • C, with other factors kept constant (e.g., with 12 h:12 h alternation of light and dark and good gas exchange). The experiments testing the effects of light and pollutant stresses were performed at a constant temperature of 25 • C. Light treatments were maintained in consistently illuminated (LL), completely dark (DD), and 12 h:12 h alternation of light and dark conditions (LD).
The effects of pollution stresses on the DRB assemblage were carried out along gradients of the pollutants. To examine the effects of antibiotics, oxytetracycline, one of the most widely used antibiotics in farming industrials, was used in this study. The stock solution (100 mg/L) of oxytetracycline (OTC, Sigma, St. Louis, MO, United States) was added, drop-by-drop, to the ciliate cultures and gently mixed, resulting in final concentrations of 1, 10, and 20 mg/L (hereafter referred to as OTC1, OTC10, and OTC20, respectively). Three types of heavy metals, lead chloride (PbCl 2 ), cadmium chloride (CdCl 2 ), and mercury chloride (HgCl 2 ; AR grade; Sinopharm, Shanghai, China) were used at final concentrations as follows: 10, 30, and 50 mg/L for PbCl 2 (hereafter referred to as Pb10, Pb30, and Pb50); 0.01, 0.05, and 0.1 mg/L for HgCl 2 (hereafter referred to as Hg0.01, Hg0.05, and Hg0.1); and 0.1, 0.5, and 1 mg/L for CdCl 2 (hereafter referred to as Cd0.1, Cd0.5, and Cd1). Since neither P. bursaria in the Cd treatments nor E. vannus in the Hg treatments survived, these treatments were not included in the subsequent analyses. The cultures without any amendments of antibiotic and heavy metals (i.e., the treatments OTC0, Pb0, Hg0, and Cd0) were taken as the controls.
For molecular assays of DRB inside the protistan cells, approximately 20-30 ciliate cells in all of these treatments were isolated during their exponential growth stages, which usually took place within 8 days of cultivation according to our pilot studies. In order to minimize the contamination from the bacteria attached to the cell surface and the cilia, protistan cells from these treatments were washed repeatedly (at least three times) by transferring to a series of sterilized seawater (for E. vannus) or autoclaved Milli-Q water (for P. bursaria) with a clean micropipette. These cleaned ciliates were then maintained in the sterilized water for two or three days, allowing the ingested bacteria to be digested as much as possible. The remaining bacteria inside these ciliate cells were considered DRB. Since the free-living bacteria in the culture medium were apparently sources for the DRB, the diversity of the former might be altered under different conditions, which was a potential mechanism modulating the diversity of the latter. To investigate this, therefore the ciliate-free culture medium was also collected for analyzing the free-living bacterial community composition and structure in each treatment.

DNA Extraction, Clone Libraries, and Sequencing
Three to five ciliate individuals with a minimum volume of water were transferred to a microfuge tube for DNA extraction using a REDExtract-N-Amp Tissue PCR Kit (Sigma, St. Louis, MO, United States), according to the modified protocol (Gong et al., 2014). Furthermore, the DNA of bacterioplankton was prepared by taking 3 µL of ciliate-free culture medium from each treatment using the same method.
The 16S rRNA gene sequences (about 1300 bp in length) obtained from cloning and sequencing were aligned using MAFFT v.7.0 with default parameters (Katoh and Standley, 2013). Potential chimeras were identified using Bellerophon (Huber et al., 2004) and Mallard 1.02 (Ashelford et al., 2006) and discarded before subsequent analysis. The cleaned sequences were clustered into operational taxonomic units (OTUs) at 97% identity using Mothur v.1.17 (Schloss et al., 2009). Taxonomic identities and systematic placements of the OTUs were assigned using the RDP Classifier (release 11.5) (Wang et al., 2007). The 16S rRNA gene sequences have been deposited in the GenBank databases with accession numbers MH555909 to MH556816.

Phylogenetic Analyses
For several sequences that assigned into unclassified Bacteroidetes and Verrucomicrobia by the RDP classifier, phylogenetic analyses were carried out using both maximum likelihood (ML) and Bayesian inference (BI) algorithms to further resolve their taxonomic ranks. Reference sequences of previously described species were retrieved from the GenBank database and aligned with the sequences in question using MAFFT. The alignment was refined manually using SeaView4 (Galtier et al., 1996). The models of nucleotide substitution were estimated using jModeltest2 (Darriba et al., 2012) based on Akaike information criterion (AIC) and Bayesian information criterion (BIC). The ML analysis was conducted using PhyML v.3.3 (Guindon et al., 2009) with a best-fit GTR + I + G evolutionary model and 1000 bootstrap replicates. The BI trees were generated using MrBayes 3.2.7 with selected parameters (rates = invgamma, number of substitution types = 6, ncat = 4) corresponding to the best model selected (GTR + I + G) (Ronquist and Huelsenbeck, 2003), implemented using Markov Chain Monte Carlo simulations and four chains running for 3,000,000 generations. The first 25% of generations were discarded as burn-in.

Profiling of Bacterial Assemblages
Variations in the assemblage structure of both the DRB associated with the protists and the bacterioplankton in the culture medium (but outside the ciliate cells) among the treatments were assessed using terminal restriction fragment length polymorphism (T-RFLP). Bacterial 16S rRNA genes were amplified as described above, except that the forward primer 63F was labeled with a fluorescent dye, 6-carboxyfluorescein, at its 5 end. The PCR products were purified with the TIAN Quick Midi Purification Kit (Tiangen, Beijing, China). The DNA yield was quantified using a NanoDrop 2000C spectrophotometer (Thermo Fisher, Wilmington, DE, United States). A 15 µL digestion mixture containing 7.5 µL of purified PCR product (each with about 20 ng of DNA) and 7.5 µL of digestion buffers with restriction endonuclease HhaI (Fermentas) were incubated at 37 • C for 1 h. Sizes of the terminal restriction fragments (T-RFs) were determined using an ABI 3730 DNA analyzer (Sangon Biotech, Shanghai, China). The baseline threshold for signal detection was set to 50 fluorescence intensity units to eliminate any background interference. Only the peaks with T-RF lengths ranging from 40 to 400 bp were included in the subsequent analysis. The relative abundance of each T-RF was calculated as the peak area ratio of the T-RF to all T-RFs detected for a given sample. Minor peaks with a relative abundance of < 1% of the total were excluded, and the remaining peaks were presumed to represent bacterial phylotypes.
The assemblage information generated from the T-RFLP profiling was analyzed using the software package PRIMER v.6 (Primer-E, Plymouth, United Kingdom). In order to test the differences in DRB assemblages, relative abundances of T-RFs were log-transformed, and a matrix of assemblage similarities was generated using the Bray-Curtis index. Nonmetric multidimensional scaling (MDS) was executed to visualize the differences in DRB assemblage structure among the treatments of each factor. To test the hypotheses that a given factor would affect the assemblage structure of the DRB and the bacterioplankton in the waters, analysis of similarity (ANOSIM) was performed.

Community Profiling of Bacterioplankton in the Culture Media
The T-RFLP profiling and ANOSIM results of global tests showed that all environmental factors investigated in this study, i.e., temperature (R ≥ 0.71, P ≤ 0.014), light (R ≥ 0.38, P ≤ 0.039), OTC (R ≥ 0.47, P = 0.001), Pb (R ≥ 0.75, P ≤ 0.002), Hg (R = 1.00, P = 0.001), and Cd (R = 0.76, P = 0.001), significantly affected the community structure of bacterioplankton in the culture media for the tested ciliate species (Table 1). From the point view of the R statistic values of global tests, the effects of the tested heavy metals on the bacterioplankton community structure seem to be stronger than those of temperature, OTC antibiotics, and light.

Community Profiling of DRB
The DRB assemblages were clearly separated into three groups by temperature in E. vannus, which was well-illustrated in the MDS plots ( Figure 1B), while the temperature-specific clusters somewhat overlapped in P. bursaria ( Figure 1A). The analysis of the assemblage structure of the DRB generally showed similar results of the bacterioplankton, namely, the effects of temperature (global test, R ≥ 0.50, P = 0.001), OTC (R ≥ 0.49, P = 0.001), and Pb (R ≥ 0.80, P = 0.001) were significant for the DRB in both species. Testing on a single species also illustrated that Hg (R = 0.73, P = 0.001) was an influential factor in the assemblage structure of DRB in P. bursaria, and Cd for that in E. vannus (R = 0.79, P = 0.001). The only exception was the light treatments in P. bursaria, in which the assemblage of DRB was not significantly different among the LL, LD, and DD treatments (R = 0.05, P = 0.362) ( Figure 1C and Table 1).
The MDS plots showed that, within the same range of a given environmental factor, the DRB assemblages were clearly separated into three groups in the marine species but overlapped somewhat in the freshwater form (Figure 1). These patterns were also supported by the results of ANOSIM: the R-value of the DRB assemblage was consistently lower than that of the bacterioplankton community in all treatments of P. bursaria (i.e., R mean = R DRB -R water = -0.22; n = 5) but higher in those of E. vannus (i.e., R mean = R DRB -R water = 0.15; n = 5; Table 1). This indicated that the DRB assemblage was less variable in P. bursaria but more variable in E. vannus as compared with their surrounding bacterial prey under varying temperature, light, OTC, Pb, Hg, and Cd concentrations. The DRB and bacterioplankton community were also significantly different in a given treatment, especially under highly stressed conditions (e.g., at 25 • C or in OTC20, see Supplementary Table S1).
Verrucomicrobia (11.3%), and Firmicutes (0.53%) were relatively minor. The proteobacterial DRB belonged to proteobacterial classes (Gammaproteobacteria, Alphaproteobacteria, and Betaproteobacteria), of which Gammaproteobacteria was the most abundant in P. bursaria and E. vannus, accounting for 37.1 and 71.9%, respectively. The DRB of Alphaproteobacteria showed similar proportions in these two ciliates (10.3 and 13.4%), whereas the relative abundance of Betaproteobacteria was much higher in the P. bursaria (13.8%) than in E. vannus (2.1%) (Figure 2). A total of 84 OTUs were obtained for the DRB of these two ciliate species, of which 72 OTUs showed high sequence identity (97-100%) with described species, and the remaining 15 OTUs had sequence similarities lower than 97% ( Table 2). The Proteobacteria OTUs were prevalent in the DRB assemblages of both ciliates, accounting for 81% of the OTUs observed (Figure 2  and Table 2). Gammaproteobacteria (41) also dominated in terms of OTU number, which was about three times the amount of Alphaproteobacteria (14) and Betaproteobacteria (13).
The most frequently occurring OTUs from DRB in the Alphaproteobacteria were affiliated with the family Rhodobacteraceae (5 OTUs, with five genera Thalassococcus, Nautella, Ruegeria, Paracoccus, and Sulfitobacter) ( Table 2). The phylotypes of the family Caulobacteraceae were represented by 4 OTUs and two genera, Caulobacter and Asticcacaulis, which were abundant in the DRB assemblages of P. bursaria under the light treatments ( Table 2). Two OTUs belonged to the genus Sphingomonas and occurred in both ciliates. The remaining three alphaproteobacterial OTUs were Methylobacterium, Phyllobacterium, and Agrobacterium, which were all associated with the P. bursaria. Nearly all of the Betaproteobacteria OTUs were recovered from P. bursaria, except for one OTU (Pelomonas sp.) detected in both two ciliate species (Table 2). Other OTUs were affiliated with the genera Aquabacterium, Curvibacter, Hydrogenophaga, Janthinobacterium, Methylophilus, Polynucleobacter, Limnobacter, Ralstonia, and Variovorax ( Table 2).
There were 10 OTUs of DRB affiliated with the phylum Bacteroidetes. These included two OTUs of Arcicella and another two of Flectobacillus, both of which are members of the family Cytophagaceae ( Table 2). Three OTUs of Flavobacteriaceae were closely related to Ochrovirga and Flavobacterium. Two OTUs were members of Cryomorphaceae: one was affiliated with Fluviicola, and the other could represent a new genus within this family (Figure 3A), exhibiting a 90% similarity with the 16S rRNA gene of Vicingus serpentipes ( Table 2). The OTU sharing an 88% sequence similarity with Crocinitomix catalasitica was confirmed to be a member of the order Flavobacteriales but remained unresolved at the family level in the ML and Bayesian trees ( Figure 3A). The remaining OTU of this phylum was closely related to Algoriphagus ( Table 2).
The actinobacterial DRB were represented by two OTUs affiliated with Micrococcus and Propionibacterium ( Table 2). Another two OTUs were related to Anaerocolumna and Peptoniphilus (phylum Firmicutes), which were only detected in the freshwater ciliate ( Table 2). A single verrucomicrobial OTU, which was associated with the marine ciliate E. vannus, shared a low sequence similarity (89%) with Ruficoccus amylovorans (family Puniceicoccaceae), likely representing a new genus ( Figure 3B and Table 2).

Warming Effect on the Assemblage Structure of DRB
Analysis of the clone libraries indicated that the DRB assemblage in P. bursaria was dominated by Pseudomonas spp., of which the relative abundance gradually increased from 50% to 65% and 76% with increasing temperature (Table 2 and Figure 4A). Flavobacterium, Peptoniphilus, and Micrococcus, though not abundant, were only detected in the warmest treatment for the freshwater ciliate. In contrast, the DRB assemblage in E. vannus was generally dominated by Pseudoalteromonas spp., of which the relative abundance decreased from 74% to 40% with temperature. The relative abundance of Vibrio sp. was minor (2.7%) during the treatment at 16 • C but became abundant (50%) at 25 • C ( Table 2  Notes: Consistently illuminated, completely dark, and 12 h:12 h alternation of light and dark conditions referred to as LL, DD and LD; oxytetracycline concentrations of 0, 1, 10, and 20 mg/L referred to as OTC0, OTC1, OTC10, and OTC20, respectively; lead concentrations of 0, 10, 30, and 50 mg/L referred to as Pb0, Pb10, Pb30, and Pb50, respectively; mercury concentrations of 0, 0.01, 0.05, and 0.1 mg/L referred to as Hg0, Hg0.01, Hg0.05, and Hg0.1, respectively; cadmium concentrations of 0, 0.1, 0.5, and 1 mg/L referred to as Cd0, Cd0.1, Cd0.5, and Cd1, respectively. and Figure 4A). Ruficoccus-like species (Verrucomicrobia) were consistently present, accounting for a small portion (3.3-13.5%) of the DRB in association with E. vannus. Nevertheless, warming induced a gradual decrease in its relative abundance.
Aestuariibacter and Pseudomonas spp. were present only in the warmest treatment of the Euplotes (Figure 4A).

Variations in the Bacterial Assemblages Under Photoperiod Conditions
E. vannus under the condition of LD hosted 9 bacterial genera, which was much more diverse than that in both LL and DD. Structurally, Pseudoalteromonas spp. and Psychrosphaera saromensis dominated under conditions of both LL and DD, whereas Aestuariibacter spp. were the most abundant (31%) in the LD treatment. Phylotypes of Aestuariibacter, Propionibacterium, Sphingomonas, Spongiibacter, and Ruficoccuslike were detected only under the LD treatment ( Figure 4B). In contrast, P. bursaria consistently hosted Nevskia sp. (38-44%), Curvibacter (25-29%), and Asticcacaulis (6-25%). Interestingly, Caulobacter and Variovorax were present but not in LD and DD, respectively ( Figure 4B).

Bacteria Co-resistant to Antibiotic and Digestion
Sphingomonas and Arcicella species occurred in the DRB of P. bursaria, accounting for 73% in the control (OTC0). However, these two taxa were not detectable in the OTC treatments. Pseudomonas in the control P. bursaria had a low abundance (19%) but a high abundance (50%-100%) after OTC applications. Also, Sulfitobacter was absent in the control and at a lower dose of OTC but represented 47% of the DRB in the OTC20 treatment ( Figure 5A). The bacteria associated with E. vannus treated with OTC included Aestuariibacter, Acinetobacter, Alcanivorax, Alteromonas, Pseudomonas, Nautella, Thalassococcus, Sphingomonas, Spongiibacter, and Ruficoccus-like Verrucomicrobia ( Figure 5A). In contrast to the non-detection of Sphingomonas in P. bursaria at high doses of OTC, this genus occurred with high relative abundance in E. vannus treated with 10 and 20 mg/L of OTC. Both the Aestuariibacter and Ruficoccuslike phylotypes were abundant at the low-dose treatment but were reduced or disappeared at the higher doses ( Figure 5A).

DRB Assemblages Changed With the Concentrations of Heavy Metals
The DRB in the control P. bursaria (Pb0) were represented by nine bacterial genera, of which only four were detected in the Pb-treated treatments, including Flectobacillus, Pseudomonas, Fluviicola, and Serratia. Flectobacillus was the most dominant, accounting for 83%, 100%, and 64% in Pb10, Pb30, and Pb50, respectively. Fluviicola and Serratia together represented 36% sequences in Pb50 (Table 2 and Figure 5B).
The DRB assemblage in E. vannus was consistently dominated by a single genus, Aestuariibacter, in both Pb and Cd treatments (Table 2 and Figures 5B,C). Nevertheless, Pelomonas (11%) occurred in Pb30 and increased to 17% in Pb50 (Figure 5B). Furthermore, in the treatments of Cd0.1 and Cd1.0, Pelomonas was also present but with minor proportions (3 and 8%, respectively). Ruficoccuslike Verrucomicrobia was minor (3.6%) in the control but abundant (∼35%) in both Cd0.1 and Cd0.5 treatments ( Table 2 and Figure 5C). FIGURE 4 | Variations in the assemblage structure of digestion-resistant bacteria (DRB) in association with two ciliates Paramecium bursaria and Euplotes vannus in different treatments. (A) A temperature gradient (16 • C, 21 • C, and 25 • C); and (B) light conditions (LL, DD, and LD represent consistently illuminated, completely dark, and 12 h:12 h alternation of light and dark conditions, respectively). Note that Pseudomonas and Pseudoalteromonas dominated the DRB assemblages in P. bursaria and E. vannus, with increasing and decreasing relative abundance with temperature, respectively. In different light environments, the composition of DRB in P. bursaria was rather similar, whereas Aestuariibacter and Alteromonas in E. vannus became more important under the condition of light-dark diel cycle. Abbreviations: Actinobacteria (Actino), Alphaproteobacteria (Alpha), Bacteroidetes (Bact), Betaproteobacteria (Beta), Gammaproteobacteria (Gamma), and Verrucomicrobia (Verruco).

DISCUSSION
By using E. vannus and P. bursaria as models, this study is the first to explicitly investigate how the composition and structure of DRB assemblages vary along various physical and chemical gradients. Overall, our molecular profiling of these two protist species indicated that experimentally manipulated warming, irradiance, antibiotic, and heavy metal stresses significantly influenced the structure of DRB assemblages. This provides evidence that the association between DRB and protists is highly dependent on environmental conditions. It was conceivable that there were consistently larger variations in the assemblage structure of DRB, relative to those in the environmental bacterioplankton across the treatments of E. vannus (Table 1). This is likely related to the fact that the variations in the environmental factors caused the shifts in the bacterioplankton community structure first (Supplementary Table S1), and then selective ingestion and differential digestion of the engulfed bacterial taxa took place next, which further screened out digestible populations, resulting in the altered assemblage structure of the bacteria (i.e., DRB) that resided in the protistan cells (Table 1 and Figures 4, 5). However,  (Hg0, Hg0.01, Hg0.05, and Hg0.1) for P. bursaria at 0 mg/L (Hg0), 0.01 mg/L (Hg0.01), 0.05 mg/L (Hg0.05), and 0.1 mg/L (Hg0.1), respectively; (C) and under cadmium treatments (Cd0, Cd0.1, Cd0.5, and Cd1) for E. vannus at 0 mg/L (Cd0), 0.1 mg/L (Cd0.1), 0.5 mg/L (Cd0.5), and 1 mg/L (Cd1), respectively. Note that Pseudomonas, Sulfitobacter, and Sphingomonas were highly antibiotic-resistant. Flectobacillus and Fluviicola were highly Pb-resistant. Hydrogenophaga was Hg-resistant, and both Aestuariibacter and Pelomonas were co-resistant to Pb and Cd. Abbreviations: Actinobacteria (Actino), Alphaproteobacteria (Alpha), Bacteroidetes (Bact), Betaproteobacteria (Beta), Gammaproteobacteria (Gamma), and Verrucomicrobia (Verruco). this reasoning cannot explain the results for P. bursaria, in which the assemblage structure of DRB across multiple environmental gradients appeared to be less variable than that of environmental bacterioplankton. The contrasting variability of DRB in P. bursaria and E. vannus may be related to their nutritional models (mixotrophy vs. heterotrophy), in which the mixotrophs have a highly flexible physiology to cope with the changing environment by shifting between phototrophic and phagotrophic lifestyles. The characteristic intracellular environment of mixotrophs thus selects a narrow spectrum of bacterial species to reside intracellularly.
The selection of DRB in all the treatments investigated in this study may be related to the direct influence of environmental factors and pollution stresses on the physiology and biochemistry of the protists, which may affect the interactions between the predators and preys. At a higher temperature (not higher than the optimal temperature), a protist usually has a high growth rate (Montagnes and Weisse, 2000), ingestion rate (Izaguirre et al., 2012), and digestion rate (Sherr et al., 1988). Compared to the light condition, protozoan ingestion rate determined by using fluorescently labeled prey was lower under the dark conditions (Izaguirre et al., 2012). Light availability is important for mixotrophs. For example, Chlorella symbionts provide photosynthetic products to their host under light conditions (Shibata et al., 2016), but they lost their ability to resist host digestion in darkness (Kodama and Fujishima, 2014). Heavy metals also affect protistan ingestion rates (Al-Rasheid and Sleigh, 1994) and physiology .
The bacterial species associated with the protistan grazers may be related to not only the life stage and physiology of ciliates , but also depend on the biological nature (e.g., digestibility) of these bacteria (Zou et al., 2020), which are more or less reflected by their phylogenetic classification. Our analyses of clone libraries revealed that the overall DRB assemblage composition in E. vannus across all treatments was dominated by Gamma-and Alphaproteobacteria, which was consistent with previous reports for a range of marine ciliate species and picoeukaryotes (Farnelid et al., 2016;Gong et al., 2016). However, this compositional trait of DRB was not applicable for the DRB in the freshwater species P. bursaria, in which Gammaproteobacteria, Bacteroidetes, and Betaproteobacteria were the most abundant three taxa. The high proportions of Bacteroidetes and Betaproteobacteria could be due to their generally high abundance in freshwater systems (Newton et al., 2011). The phylum Bacteroidetes is well known for its high abundance in the gut microbiota (Mahowald et al., 2009) and for having a special protein secretion system (Type IX Secretion System, T9SS) that ensures cell survival and fitness in response to its habitat (Lasica et al., 2017). It is likely that habitat (freshwatermarine) distinctness in higher taxonomic composition of DRB exists, but it is a notion that needs to be further investigated for more freshwater protist species. At lower taxonomic ranks, the DRB assemblage composition and structure along different environmental gradients and pollutant stresses showed many interesting characteristics, many of which are reported for the first time in this study, and the underlying mechanisms are discussed below.

Warming Selected Specific Bacterial Species in Association With Protists
It is interesting to observe that Pseudomonas spp., as a dominant group in P. bursaria, had higher relative abundances at higher temperatures ( Figure 4A). It was possible that the endosymbiotic green algae Chlorella sp. became more carbon-rich at higher temperatures (Serra-Maia et al., 2016), which could supply more organic matter to the symbiotic Pseudomonas species (Yao et al., 2019). An opposite trend was observed for the dominant DRB group in the algae-free ciliate E. vannus, in which the relative abundance of Pseudoalteromonas species decreased progressively with the increase of temperature ( Figure 4A). The cellular chemical composition of E. vannus could be altered with temperature, with decreasing carbon to nitrogen ratio or phosphorus ratio at higher temperatures, according to the ecological stoichiometry theory (Woods et al., 2003;Fu and Gong, 2017). This temperature-driven stoichiometric shift in the host could select against specific heterotrophic DRB populations, such as Pseudoalteromonas spp. The lower relative abundance of Pseudoalteromonas occurred at a higher temperature (25 • C), which could be attributed to the mechanism whereby Pseudoalteromonas kills other bacterial competitors, such as Vibrio, by secreting pseudoalterin, a thermolabile protease (Richards et al., 2017;Tang et al., 2020), which might lower the killing activity at a higher temperature, resulting in a higher proportion of co-occurring Vibrio in the DRB assemblage. Also, a higher temperature might enrich the type III secretion system (T3SS) of Vibrio in association with protists (Matz et al., 2011).

Light-Driven DRB Assemblage in the Heterotrophic Protist, but Light-Independent DRB in the Mixotroph
The DRB assemblage significantly responded to the light conditions, which occurred only in the heterotrophic E. vannus but not in the mixotrophic P. bursaria, indicating that the mixotroph could select specific bacterial populations. Interestingly, the DRB detected in the P. bursaria, i.e., the gammaproteobacterial Nevskia and the alphaproteobacterial Caulobacter and Asticcacaulis, share a common morphological character in having a stalk or holdfast at certain life stages, with which cells are able to attach onto a water-air interface or substrate surface (Stürmeyer et al., 1998;Purcell et al., 2007;Persat and Gitai, 2014). Typical species of the former two genera have also been demonstrated to be light-sensitive. For example, Nevskia ramosa often inhabit at the air-water interface and is capable of DNA reparation against UV damage (Stürmeyer et al., 1998). Visible light can act to regulate a sensory module in Caulobacter crescentus in order to increase the biochemical activity and cellular signaling for cell-surface and cell-cell attachment (Purcell et al., 2007). Asticcacaulis was detected in the microbial consortium of the oil-rich green alga Botryococcus braunii (Sambles et al., 2017). In the constant light treatment, continuous uptake of nutrients by the autotrophs, which included those in medium and the Chlorella symbionts inside the P. bursaria, probably lead to nutrient limitation to heterotrophic bacteria in the medium. Attachments to, or residue of these bacteria inside the ciliate cell may facilitate themselves to immediately "capture" the nutrients regenerated by the protist. In constant darkness, P. bursaria might sustain by bacterivory, and the bacteria-protist association might reduce bacterial mortality. Therefore, our study based on light manipulation and single cell analysis of protists provides evidence that mixotrophic protists are previously unrecognized hosts and micro-niches of these attaching bacteria in natural environments.
The Curvibacter phylotype occurred as DRB of the P. bursaria in all-light, all-dark, and diel-cycle treatments. Cultivation in constant darkness or at night might lower the concentration of dissolved oxygen, which provides an ecological advantage to microaerobic bacteria, such as Curvibacter sp. (Ding and Yokota, 2010). It has been demonstrated that members of this genus are able to inhibit fungal infection by interacting with other commensal bacterial species in the cnidarian Hydra (Fraune et al., 2015). In fact, at least filament fungi were recognizable via the naked eye at later periods of the cultures supplied with rice grains, indicating a risk of fungal infection of the P. bursaria. In this sense, it is possible for these bacterial consortiums to function in anti-fungal activity. The Variovorax phylotypes were only detected in both LL and LD treatments of the P. bursaria. This adaptive ability in such contrasting light environments is consistent with a previous study, which showed a combination of autotrophic and heterotrophic features in the genome of Variovorax paradoxus (Han et al., 2011).
The DRB in the heterotrophic ciliate E. vannus were photoperiod-dependent. Alteromonas sp. was dominant in the treatment of the diel light:dark cycles, suggesting a role of light (or darkness) in selecting these bacteria in association with the protistan cell. A similar mechanism was proposed by a previous study (Biller et al., 2018), which stated that Alteromonas may provide energy or organic compounds to cyanobacteria Prochlorococcus under the stress of darkness. The Pseudoalteromonas phylotypes were consistently present in the DRB assemblages under the conditions having a photic period, probably contributing to anti-fungal infection by more actively producing antifungal polyketide alteramides during the dark period (Moree et al., 2014). Thus, the co-occurrence of these two bacterial taxa may be of benefit to the host for a balance between rapid growth and low pathogenic infection under a natural diel light-dark condition (Liu et al., 2019). Gomiero and Viarengo (2014) showed that the antibiotic OTC (3.2 µM to 32 mM) treatments induced changes in the growth, survival, endocytosis rate, and lysosomal membrane stability of the marine ciliate Euplotes crassus. Such physiological shifts might have also occurred in these two species we investigated. Apart from that, we found the dominant DRB in the OTCtreated marine ciliate E. vannus were affiliated with the members of the genus Sphingomonas and the family Alteromonadaceae (Alteromonas and Aestuariibacter) (Figure 5A), which is in line with previous findings that both Sphingomonas and Alteromonas species were resistant to antibiotics, such as OTC (Miranda and Zemelman, 2002;Dang et al., 2007). However, the relative abundance of the phylotype of the family Puniceicoccaceae (phylum Verrucomicrobia) decreased in OTC10 and was nondetectable in OTC20, indicating it is sensitive to OTC.

Evidence for Prevalent Antibiotic-Resistant Bacteria Sheltered in Protistan Cells
Both Pseudomonas and Sulfitobacter phylotypes prevailed in the DRB assemblage of mixotrophic P. bursaria reared in media with an antibiotic concentration up to 20 mg/L, indicating high antibiotic resistance of these bacteria and that the intracellular environment may provide a favorable habitat for their survival. Indeed, it has been demonstrated that many Pseudomonas species could degrade antibiotics via efflux pumps (e.g., Jiang et al., 2014), and consortiums of these bacterial species with algae could enhance the efficiency of degrading antibiotics (Leng et al., 2020). Similarly, genes involved in aromatic compound catabolism and a type IV secretion system are present in the genome of a strain of Sulfitobacter (Ankrah et al., 2014), and Sulfitobacter in the diatom phycosphere could supply indole acetic acid to the algae in exchange for organosulfur compounds (Amin et al., 2015). Furthermore, Sulfitobacter strains isolated from marine hydrothermal vent fields were found to be antibiotic-resistant (Farias et al., 2015). All these suggest that the mixotrophic protist bearing endosymbiotic microalgae (e.g., Chlorella sp.) could play a role in accommodating these antibiotic-resistant bacteria, which represent a previously unrecognized scenario in the antibiotic resistance of environmental bacterial populations.

Bacterial Co-resistance to Heavy Metals and to Protistan Digestions
It is well known that short-rod cells of Flectobacillus form long filaments and chains of several cells, a morphology-based strategy to protect themselves from being ingested by nanoflagellates predators (Corno and Jürgens, 2006). We detected these bacteria as a member of DRB in P. bursaria, suggesting an effect of predator size in the prey-predator interactions, i.e., enlarged bacterial cell size may be effective to avoid being engulfed by small protists but not necessarily for large-celled protists, such as ciliates with wide or highly contractile cytostomes. Furthermore, the Flectobacillus phylotypes became dominant in the Pb-treated cultures, demonstrating they are co-resistant to Pb and protistan digestion. However, Flectobacillus spp. were found to be a dominant bacterial group in the riverbed sediments spiked with high concentrations of metals, including Cd and Cu, but not in the Pb and Cr treatments (Du et al., 2018). In the Hg-treated P. bursaria culture, we also found that Flectobacillus phylotypes only accounted for a minor portion of DRB. This suggests a strain-level differentiation in the heavy metal resistance of Flectobacillus.
Acinetobacter junii is able to form biofilms on surfaces to resistant to mercury (Sarkar and Chakraborty, 2008), which is in line with our observation of its survival in P. bursaria under Hg stresses. In addition, our results point to the non-digestible nature of this bacterial species in ciliated protozoa, highlighting its ecological success in natural environments. Furthermore, Hydrogenophaga was also abundant in P. bursaria treated with a high dose of Hg. It remains to be investigated how this hydrogenoxidizing species interacts with the host and endosymbiotic green alga within the P. bursaria under Hg stress, though it is known that As (III) could be oxidized by a strain of this genus under aerobic conditions (Terry et al., 2015).
Not much is known about the ecology and microbial association of Aestuariibacter, except for our recent report of this genus as a DRB in four marine and one freshwater species of ciliates (Gong et al., 2016). In the present study, we once again found this genus as a dominant group of DRB in E. vannus across all Pb 2+ and Cd 2+ treatments and under low-dose antibiotic treatment, indicating Aestuariibacter is a previously unrecognized superstar in co-resistance to heavy metals, antibiotics, and protistan digestion. This is probably due to their resident habitats of high-pollution stresses, such as coastal seawater and tidal flats (Yi et al., 2004;Wang et al., 2010). Pelomonas was present in the DRB of E. vannus treated with high doses of Pb 2+ and Cd 2+ but absent in the controls, indicating the metal tolerance of this group. Phylotypes of this genus were also detected in Cd-and Cu-treated activated sludge using DGGE and sequencing (Bhat et al., 2020). Strains of Verrucomicrobia have been isolated as gut symbionts from sea cucumber, marine sponges, clamworm (Wertz et al., 2012), and the marine ciliate Euplotidium (Petroni et al., 2000). Our study showed the association between a possible new genus of the family Puniceicoccaceae with the ciliate E. vannus and for the first time their differential resistance to metals Pb and Cd.

CONCLUDING REMARKS
As major players in the microbial loop and biogeochemical cycle within ecosystems, both bacteria and protists are ecologically linked not only through prey-predator food chains but also via collaboration in facing multiple stresses. Our study reveals for the first time that the assemblage composition and structure of ingested but inedible bacteria in two protists vary significantly in response to stresses of temperature, light, antibiotics, and heavy metals, indicating that the environment plays an important role in selecting specific bacterial populations in protist-bacteria associations. Our results for these two ciliate species also provide an indication that trophic mode of protists (mixotrophs vs. heterotrophs) could be a factor influencing the bacteriaprotist interactions, a notion warrants further investigation considering their increasingly recognized importance in marine food webs (Stoecker et al., 2017). Furthermore, bacterial antibiotic resistance and metal resistance have been extensively studied in the past several decades. The ecological interactions between these resistant bacteria and protists remain poorly understudied (Nguyen et al., 2020). Our findings in this study provide evidence that there are diverse bacterial stains, which are not only resistant to antibiotics and/or heavy metals in natural environments but have also survived protistan digestion, suggesting some bacterial species/strains could use protistan cells as shelters in facing of pollution stresses.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repositories and accession numbers can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
SZ did the data curation-equal, formal analysis-equal, investigation-equal, methodology-equal, resources-equal, software-equal, validation-equal, visualization-equal, and wrote the original draft-lead. QZ performed the resources-supporting and validation-supporting. XZ and CD performed the validationsupporting. JG performed the conceptualization-lead, data curation-supporting, funding acquisition-lead, methodologyequal, project administration-lead, resources-equal, supervisionequal, validation-equal, and wrote, reviewed and edited the manuscript-lead. All authors contributed to the article and approved the submitted version.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2020.00659/full#supplementary-material FIGURE S1 | A general scheme of research methodology.