ORIGINAL RESEARCH published: 28 April 2022 doi: 10.3389/fmicb.2022.825167 Recent Differentiation of Aquatic Bacterial Communities in a Hydrological System in the Cuatro Ciénegas Basin, After a Natural Perturbation Manuel García-Ulloa II 1,2 , Valeria Souza 1,3* , Diego A. Esquivel-Hernández 1 , Jazmín Sánchez-Pérez 1 , Laura Espinosa-Asuar 1 , Mariette Viladomat 1 , Montserrat Marroquín-Rodríguez 4 , Marisol Navarro-Miranda 1 , Jair Ruiz-Padilla 1 , Camila Monroy-Guzmán 1 , David Madrigal-Trejo 1 , Manuel Rosas-Barrera 1 , Mirna Vázquez-Rosas-Landa 1 and Luis E. Eguiarte 1 1 Departamento de Ecología Evolutiva, Instituto de Ecología, Universidad Nacional Autónoma de México, Ciudad de México, Mexico, 2 Departament de Genètica i Microbiologia, Facultat de Biociències, Universitat Autònoma de Barcelona, Barcelona, Spain, 3 Centro de Estudios del Cuaternario de Fuego-Patagonia y Antártica (CEQUA), Punta Arenas, Chile, 4 Facultad de Medicina, Licenciatura en Investigación Biomédica Básica, Universidad Nacional Autónoma de México, Ciudad de México, Mexico Edited by: Céline Lavergne, Universidad de Playa Ancha, Chile Pozas Rojas is a hydrological system comprising nine isolated shallow ponds and Reviewed by: a deep lagoon, which were temporally merged in 2010 by increased rainfall due to Janet Siefert, a tropical cyclone. In this work, we assess which components, biotic interactions, Rice University, United States Polette Aguilar-Muñoz, or environment filtering effects, drive the assembly of microbial communities after a Universidad de Playa Ancha, Chile natural perturbation. Arsenic, pH, and temperature are among the most significant *Correspondence: environmental variables between each pond, clustering the samples in two main groups, Valeria Souza souza@unam.mx whereas microbial composition is diverse and unique to each site, with no core at the operational taxonomic unit level and only 150 core genera when studied at the genus Specialty section: level. Los Hundidos lagoon has the most differentiated community, which is highly similar This article was submitted to Systems Microbiology, to the epipelagic Mediterranean Sea communities. On the other hand, the shallow ponds a section of the journal at the Pozas Rojas system resemble more to epicontinental hydrological systems, such Frontiers in Microbiology as some cold rivers of the world and the phreatic mantle from Iowa. Overall, despite Received: 30 November 2021 Accepted: 25 February 2022 being a sole of water body 2 years prior to the sampling, interspecific interactions, Published: 28 April 2022 rather than environmental selection, seem to play a more important role in Pozas Rojas, Citation: bolstered by founder effects on each poza and subsequent isolation of each water body. García-Ulloa M II, Souza V, Esquivel-Hernández DA, Keywords: Eltonian niche, Grinnellian niche, founder effect, environmental filtering, compositional differentiation, Sánchez-Pérez J, Espinosa-Asuar L, ecological perturbation Viladomat M, Marroquín-Rodríguez M, Navarro-Miranda M, Ruiz-Padilla J, Monroy-Guzmán C, Madrigal-Trejo D, INTRODUCTION Rosas-Barrera M, Vázquez-Rosas-Landa M and Early in the 20th century, two main frameworks for the study of ecological niches have been Eguiarte LE (2022) Recent proposed: the Grinnellian and Eltonian niche concepts. On one hand, the Grinnellian niche Differentiation of Aquatic Bacterial concept (Grinnell, 1904) was meant as a description of the complete range of conditions Communities in a Hydrological System in the Cuatro Ciénegas Basin, and resources where any given organism could live and reproduce (see review in Soberón, After a Natural Perturbation. 2007). On the other hand, the Eltonian (Elton, 1927) niche concept is explained by the Front. Microbiol. 13:825167. interactions among species, with the trophic web as central for the establishment of local doi: 10.3389/fmicb.2022.825167 community assembly (Soberón, 2007). Although there are multiple studies on the interplay of the Frontiers in Microbiology | www.frontiersin.org 1 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities Grinnellian and Eltonian niches in natural animal and plant recently for Vibrio (Vázquez-Rosas-Landa et al., 2017, 2020). It populations (Soberón, 2007; Larson et al., 2010; Rosado et al., is interesting to note that earlier metagenomics of Pozas Rojas 2016; Junker et al., 2019), their role in succession processes and observed a strong dominance and diversity of Pseudomonas, a perturbation scenarios is more difficult to dissect, and questions lineage that probably benefited from the environmental filtering regarding the relative importance of local environmental due to the extreme unbalanced stoichiometry of the site (C:N:P conditions and interactions among members of the community 15,820:157:1; Bonilla-Rosso et al., 2012; Peimbert et al., 2012). remain to be answered, particularly in in situ microbial setups. However, the extreme oligotrophy of Pozas Rojas was perturbed Succession process in aquatic microbial communities is by Hurricane Alex in April 2010. This natural perturbation seldom studied in situ given the complexities of natural brought not only water that temporarily connected all pozas and communities and a large amount of confounding factors, the lagoon into one body of water, but also debris from most of including community response to environmental variables the east side of the basin to the region named Los Hundidos due to and stressors, stochastic processes, migration, and founders’ its karstic nature and slightly lower elevation than the rest of the effect (Pajares et al., 2013, 2015; Pascual-García and Bell, basin. As time passed, the superficial excess of water refilled the 2020). However, aquatic mesocosm experiments have shown deep aquifer, leaving in 2012 the pozas in the same geographical that environmental variables such as temperature and UV position they were before the perturbation. At the moment of light determine the response of the microbial community to the sampling in 2013, the site was no longer oligotrophic (C:N:P perturbation, as the Grinnellian niche concept would predict. 350:9:1) allowing the more copiotroph Vibrio to thrive (Vázquez- Other studies have shown the influence of the Eltonian niche. Rosas-Landa et al., 2020). Also, in a recent sampling of Pozas For instance, the establishment of bacterial communities in beech Rojas viral diversity, compared with Churince, the pozas showed tree holes seems to be governed by functional redundancy and a much larger diversity and differentiation in the water column founders’ effect (Pascual-García and Bell, 2020). In contrast, between water bodies in March 2014 (Taboada et al., 2018), a year communities from oligotrophic lakes in the Pyrenees have shown after the present study. to be influenced by the interactions between environment and the Herein we studied the aquatic microbial diversity of Pozas community as well as migration processes, being environmental Rojas, 3 years after Hurricane Alex merged all the pozas into filtering a strong force (Ortiz-Álvarez et al., 2020). Therefore, one lake and then got reseparated. In this work, we want both the Grinnellian and Eltonian niches are relevant for the to address the issue of differentiation between pozas in a establishment of the lake microbial communities. very short time, testing the hypothesis of Grinnellian versus Another fundamental issue is what happens to microbial Eltonian niche as drivers of the community structure. As a communities after natural perturbations, that is, when the null hypothesis, the flooding erased the structure between pozas Grinnellian niche changes abruptly. In a 4-year study of halite within Pozas Rojas, and at that time, they showed no significant communities from the Atacama Desert, after an unusual rain, differentiation. However, there are two alternative hypotheses: (1) the community shifted to a completely new unstable state the Grinnellian niche drives the pozas differentiation; therefore, and then recovered through stochastic recolonization once environmental variables explain most of Pozas Rojas diversity; (2) the water evaporated (Uritskiy et al., 2019). Opposite to the each poza will be unique within Pozas Rojas due to the Eltonian Atacama example, in the Churince system from Cuatro Ciénegas niche where interactions rule community assembly. Moreover, Basin (CCB), water became scarce, and salinity progressively we scale such geographical differentiation from the local to a rose due to the overexploitation of the aquifer. This resulted more global scale by comparing the Pozas Rojas composition in demographic fluctuations of the Pseudomonas community, with other epicontinental and marine water bodies: the phreatic producing a sudden bloom of Pseudomonas aeruginosa (García- mantle from Iowa, a groundwater system; the Mediterranean Sea, Ulloa et al., 2019) and an event of evolutionary rescue of a marine system; several cold rivers from the world, which are Pseudomonas otitidis (García-Ulloa et al., 2021). However, as the instances of lotic systems; the Churince system from the CCB; perturbation became stronger and disrupted the community, the and the well-studied Lake St. Clair in Michigan, a lentic system entire genus disappeared and a never-before-seen environmental that is part of the Great Lakes Basin. We expected sites within strain of Vibrio cholerae appeared close to the extinction of the Pozas Rojas to be more similar among each other than other CCB lagoon. The interplay between Grinnellian and Eltonian niches sites, such as Churince. Following this logic, CCB’s systems would seems to explain this outcome; as the Grinnellian niche radically be more similar between them than between different sites of the changed, the Eltonian niche eventually became altered, and the world, as Pozas Rojas and Churince may share a core of CCB original community collapsed. Unfortunately, as in the case of microbiota given their shared deep aquifer (Wolaver et al., 2013). P. otitidis and possibly most of the Pseudomonas genus, the enormous diversity of Churince [5,167 operational taxonomic units (OTUs) in a square km (Souza et al., 2018)] is now mostly MATERIALS AND METHODS lost due to the water depletion. Within the studied sites at CCB is Pozas Rojas, located Sample Collection, DNA Purification, and on the east side of the valley and composed of multiple Sequencing shallow pozas and a deeper lagoon. This place was studied in Pozas Rojas system, within the CCB, is located at the geographic parallel to Churince as a comparative site for cultivated Bacillus, coordinates 26◦ 520 N, 102◦ 10 W. Sampling site overview and Exiguobacterium, and Pseudomonas (Avitia et al., 2014) and more geographic coordinates for each poza/lagoon are provided in Frontiers in Microbiology | www.frontiersin.org 2 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities the Supplementary Figure 1 and Supplementary Table 1. The 360◦ C. Soil N was determined by the macro-Kjeldahl method pozas (designated S02–S09 in the present study) are small and (Bradstreet, 1954), whereas P was determined by the molybdate shallow (no more than 20 cm in the winter); therefore, wind colorimetric method following ascorbic acid reduction (Dick and is expected to evenly mix the water column, except at the Los Tabatabai, 1977). The N and P forms analyzed were determined Hundidos (designated LH1) lagoon, which is roughly 10 m deep, colorimetrically in a Bran–Luebbe Autoanalyzer 3 (Norderstedt, and stratification is presumed. Each water body was sampled with Germany). Several replicates per measurement were used in the sterilized water bottles, collecting 6 L of superficial water and statistical analysis. keeping them in a cold icebox. Simultaneously, environmental variables were measured with a Hydrolab MS5 multiprobe. Three Microbial Diversity liters from each water body was taken for mineral and nutrient A total of 223,937 reads were obtained with MiSeq Illumina determination (CIECO and Instituto de Geofísica, UNAM). 2 × 300 using 357F and 939R primers. Reads passing the QC Microbes were collected by filtering in situ all the sampled water filters (minimum read quality = Q20) were used to reconstruct onto sterile GF/F filters (0.2-mm nominal pore size; Whatman, the original amplicon region (450–490-bp length) by overlapping Piscataway, NJ, United States) using a Millipore filtering device. them with Flash v1.2.7 software (Magoč and Salzberg, 2011). All GF/F filters were used in order to accumulate enough biomass non-overlapping sequences were discarded. A total of 188,190 for DNA extraction. sequences were used for taxonomic annotation with the Parallel- DNA was extracted from water column samples using the META v2.4.1 pipeline (Su et al., 2014) against the Metaxa2 MOBIO PowerWater DNA Isolation kit (MoBio Laboratories, database v2.1.1 (Bengtsson-Palme et al., 2015; Escobar-Zepeda Carlsbad, CA, United States), with one modification: the volume et al., 2018). Stacked bar plots at different taxonomic levels were of PW1 solution was increased to 1.5 mL due to the high generated with ggplot2 (Wickham, 2016). The OTU abundance absorbency of the GF/F filters, obtaining several filters per poza. table was used to calculate the Good’s coverage (rarefaction) and α Therefore, a composed sample of several independent DNA diversity indexes, that is, Observed, Chao 1, Fisher, Simpson, Inv isolation steps from each site was pooled into one sample per Simpson, and Shannon, with the R Phyloseq library (McMurdie poza. The DNA was used to amplify the V4–V6 16S rRNA and Holmes, 2013). The OTU table was normalized using the gene variable regions using the primers 357F and CD[R] (Rudi metagenomeSeq method (Paulson et al., 2013), and the β diversity et al., 1997; Turner et al., 1999). Amplicon libraries (450–490- distance matrix was calculated using Bray–Curtis dissimilarity. bp length) were constructed as reported in the 16S Metagenomic Absence/presence plots were generated using the Upset R library Sequencing Library Preparation protocol from Illumina and (Lex et al., 2014) from the abundance OTU table and taxonomic sequenced on the Illumina MiSeq platform with a paired-end annotation table at the genus level. read configuration of 150 cycles at CINVESTAV-LANGEBIO, Irapuato, Mexico. Statistical Analyses Statistical analyses were conducted with R.1 Using the Measurement of Physicochemical physicochemical and environmental measurements, the optimal Parameters number of clusters was assessed, from k = 2 to k = 8, with cluster v2.1.0, and plotted with factoextra v1.0.7 R packages, respectively To determine the ionic concentration of the samples, sediment (Kaufman and Rousseeuw, 2009). Pearson correlation coefficients samples were incubated with distilled water for 19 h at were calculated to assess the associations between microbial α 25◦ C under continuous shaking. This procedure allows diversity and environmental measurements and the relationships the mobilization of the available ions within the sediment. across various environmental variables. p < 0.05 was considered Liquid samples from sediment and water were filtered with to be statistically significant. 0.22 nitrocellulose membranes. Concentrations of CO3 2− and Non-correlated variables were selected based on a correlation HCO3 − anions were determined by HCl 0.01 N titration. Cl− and plot calculated with the R package corrplot v0.84 (Wei and Simko, SO4 2− concentrations were determined by high-performance 2017), resulting in eight independent environmental variables: liquid chromatography. A 432-conductivity detector (Waters) temperature, C:N, calcium (Ca), iron (Fe), conductivity, TP, and an IC-Pack HR-Waters column were used. The mobile phase pH, and arsenic (As). However, Fe was excluded from the was a mixture of sodium borate (1.3 M) and acetonitrile (12%) at analysis because of its slight presence in only two of the pH 8.5. Ca2+ and Mg2+ concentrations were obtained by atomic samples. Euclidean dissimilarity matrices of individual and every absorption spectroscopy using a Perkin Elmer 3110 equipment. possible group of non-correlated environmental variables and a Na+ and K+ were determined by flamometry using a Corning Harvestine dissimilarity matrix of the coordinates from Pozas 400 device. For nutrient quantification, sediment samples were Rojas samples were calculated. Mantel tests were performed dried, and water samples were filtered through a Millipore 0.42- with both Spearman and Pearson correlation methods (9999 µm filter. Total carbon (TC) and inorganic carbon (IC) were permutations) using the R package vegan v2.5.6 (Dixon, 2003). determined by combustion and colorimetric detection using a In addition, a canonical correspondence analysis (CCA) of TC analyzer (UIC model CM5012, Chicago, IL, United States). the environmental and physicochemical measurements was Total organic carbon was calculated as the difference between TC and IC. For total N and total P (TP) determination, samples were acid digested with H2 SO4 , H2 O2 , K2 SO4 , and CuSO4 at 1 http://cran.r-project.org Frontiers in Microbiology | www.frontiersin.org 3 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities performed using the cca function from the R package vegan of the following chemical species: CO3 2− , HCO3 − , v2.5-6. Only variables with p < 0.05 are shown to reduce noise. SO4 2− , Cl− , Na+ , K+ , Ca2+ , Mg2+ , Cd, Pb, and As (Supplementary Table 2). Pozas were separated by such Data Acquisition of Worldwide Water environmental variables into two groups: LH1, S04, S05, Bodies and Taxonomic Annotation S07, S08, and S09 in one group and S02, S03, and S06 in Pozas Rojas raw reads were uploaded to MG-RAST and another (Figure 1). are publicly available on the site2 with the project accession Univariate Mantel tests show nine significant correlations of number mgp94066. 16S tags from water samples collected geographic location with temperature and As [As (Spearman): globally were used as references. Accession numbers of the r = 0.615, p = 0.0164; temperature (Spearman): r = 0.366, reference sites are as follows: mgp82906 (phreatic zone, Iowa), p = 0.038; As (Pearson): r = 0.818, p = 0.0071; temperature mgp10029 (eastern Mediterranean Sea), mgp89153 (Lake St. (Pearson): r = 0.467, p = 0.013], and multivariable Mantel points Clair, Michigan), mgp83892 (cold rivers, worldwide), and to temperature and pH [temperature + pH (Spearman): r = 0.404, mgp94070 (Churince, Cuatro Cienegas, Mexico). All 16S tags p = 0.037]. In addition, CCA gives significant results for pH and used in this study were processed with the MG-RAST pipeline conductivity (Table 1). (Keegan et al., 2016). Rarefaction curves of the species count were performed for all geographical sites, and taxonomic annotation was done using the Ribosomal Database Project database as reference. An abundance table was constructed with the following minimum values: e = 10e−5, percentage of identity = 60%, length = 15 bp, abundance = 1. Comparative Analysis With Samples From Around the World Bray–Curtis dissimilarity matrices of relative abundances and non-metric multidimensional scaling (NMDS) analyses of all taxonomic categories, from phylum to genus, for all sites around the world were calculated using the R package vegan v2.5.6. Homogeneity of multivariate dispersion at each taxonomic category was assessed with an analysis of variance (ANOVA) of the β dispersions of the Bray– Curtis dissimilarity matrices previously calculated. Validity of the NMDS analyses was assessed with the stress scores. A permutational multivariate ANOVA (1,000 permutations) using dissimilarity matrices and an analysis of similarity (ANOSIM) were carried out to measure significant differences between sites based on their taxonomic abundances with the FIGURE 1 | Dendrogram of environmental variables. Hierarchical K-means clustering algorithm of the R package cluster v2.1.0 was used. Two main adonis and anosim functions from vegan v2.5.6, respectively. groups resulted from the clustering of environmental variables, shown in blue In addition, a stacked bar plot of the relative abundances and yellow. of phyla from all sites was done with the R package ggplot2. Rarefaction curves for each sample are shown in Supplementary Figure 4. TABLE 1 | Canonical correspondence analysis of environmental variables. Data χ2 F stats Pr (>F) Significant variables RESULTS All Temperature, pH, 1.70192 1.0343 0.188 COT As, pH, and Temperature Are the Most conductivity, C:N, (p = 0.034) C:P, N:P, COT, NT, Important Environmental Variables in PT, CO3, HCO3, Pozas Rojas SO4, Cl, Na, K, Ca, Nutrient analysis showed that the proportion of C:N:P was, on Mg, Cd, Pb, As average, 350:9:1 for water and 258:21:1 for sediment. These Environmental Temperature, pH, 1.10427 2.2101 0.013 pH (p = 0.001) conductivity Conductivity nutrients (C, N, and P) were not significantly different (p = 0.007) between sampling points for water samples (Vázquez- Stoichiometry C:N, C:P, N:P 0.97696 1.696 0.088 C:P (p = 0.004) Rosas-Landa et al., 2020). Environmental variables include Minerals COT, NT, PT, CO3, 1.70192 1.0343 0.216 CO3 (p = 0.05) pH, temperature, conductivity, and the concentration HCO3, SO4, Cl, Na, K, Ca, Mg, Cd, 2 mg-rast.org Pb, As Frontiers in Microbiology | www.frontiersin.org 4 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities despite that the somewhat closer pozas (S02, S05, and S03) are more similar among each other in terms of community composition (Figure 2). Despite its recent recolonization, the diversity and equitativity within each site is very high as reflected by the high Shannon and Simpson indexes, with the notable exception of S02 and the less well-sampled site S07 (Table 2). There Is No Core Community in Pozas Rojas The divergent community composition of LH (Figure 2) is also evident in the unique versus shared OTUs between Pozas Rojas sites (Supplementary Figure 3), with LH having 1334 unique OTUs, whereas the next site in line, S08, has only FIGURE 2 | Dendrogram of Bray–Curtis distance between samples. 403. Surprisingly, there are no core OTUs that are shared among all sites despite their proximity and the fact that they were a connected water body 2 years prior to this Each Poza Possesses a Distinctly High sampling (Supplementary Figure 3). In order to understand such Microbial Diversity uniqueness, we assigned genera to the most abundant OTUs From a total of 188,190 reads, 5,268 OTUs were obtained at 97% (>0.1%). With this approach, we were able to identify that similarity. Rarefaction curves show that even if the number of 150 genera are shared between the different pozas (Figure 3). reads per poza is not the same, whereas LH is very well sampled On the most divergent site, LH, the marine Cyanobacteria along with S02, S08, and S09, the other pozas are subsampled Synechococcus is the most abundant, whereas it is rare in the (Supplementary Figure 2). However, microbial diversity is not shallower pozas (Figure 4). Meanwhile, most shallow pozas have related to the sequence coverage or to environmental variables. diverse proportions of the Actinobacteria Candidatus Rhodoluna, Community differentiation shows an isolated site (LH1) and two as well as a high abundance of Rhodobacteraceae along with more related clusters (Figure 2 versus 1). The Los Hundidos site fewer Cyanobacteria. has clearly a unique microbial composition. This is confirmed by Bray–Curtis β diversity that is, in general, also very high, showing Migration Is Rare Between Los Hundidos that each site is diverse and unique (Figures 2, 3 and Table 2). and the Pozas Furthermore, a Mantel test shows that this differentiation is While at higher taxonomic levels clustering of sites is generally unrelated to geographic distance (p = 0.2274, R2 = 0.03291), not clear, at the genus level, all Pozas Rojas samples group FIGURE 3 | Upset plot of shared genera between pozas. Connected points represent the different combinations of sample intersections. First vertical frequency bar displays the common genera between all ponds. The total number of genera in each pond is shown by the horizontal frequency bars. Frontiers in Microbiology | www.frontiersin.org 5 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities TABLE 2 | Microbial α diversity in the water of the Pozas Rojas hydrological DISCUSSION system after a natural perturbation. Sample Observed Chao1 Shannon ACE Simpson InvSimp Fisher All the data show that there is a strong community structure and differentiation between the individual pozas, particularly LH1 1,851 1,880.57 4.12 17.7 0.92 12.24 365.59 with the deeper and most stable site, Los Hundidos lagoon. S02 1,143 1,449.02 2.85 22.95 0.71 3.5 235.81 As we stated previously, there are two possibilities: (1) the S03 415 893.85 4.54 19.91 0.96 27.5 161.68 Grinnellian environmental variables explain most of Pozas S04 736 1,220.71 5.06 22.55 0.98 49.85 255.53 Rojas diversity and structure; and (2) each poza is unique S05 1,271 1,734.97 4.92 25.37 0.97 29.97 365.65 because of the Eltonian niche determining community assembly. S06 905 1,323.49 4.94 21.68 0.97 32.46 271.47 In this study, we observed a significant differentiation in S07 497 664.2 3.04 16.63 0.72 3.55 137.69 community structure between pozas given by environmental S08 1,841 2,358.46 5.43 28.8 0.98 42.28 501.55 variables such as temperature and pH. This is to be expected, S09 1,123 1,448.9 4.28 23.36 0.95 20.21 290.76 as these environmental variables are the most important variables in microbial community assembly in general (Hollister et al., 2010; Ren et al., 2015; Yan et al., 2017). However, together and very closely to Churince samples, except for LH, nutrients and minerals, with the exception of As that is which arranges itself as a part of the epipelagic samples of more abundant in pozas S02, S04, and S07, appear to be the Mediterranean Sea (Figure 5). Moreover, a clear difference less relevant variables for community assembly in Pozas can be appreciated on the taxonomic profiles of LH1 and Rojas. What is interesting is that LH, environmentally, is the rest of the pozas along similarities with the epipelagic part of the large inner cluster (Figure 1); this could be Mediterranean Sea samples at the phylum level, mainly a due to the drainage of the hurricane water toward it. very high relative abundance of Cyanobacteria followed by Nevertheless, such clustering is not reflected in the microbial Bacteroidetes (Figure 6). This suggests not only that the water community structure, being LH the most divergent site. source of LH is different from the shallow pozas, but also that Therefore, although there is an influence of Grinnellian superficial migration from that lagoon to the pozas around niche on community assembly, Eltonian processes seem to the lagoon is rare. play a more important role in Pozas Rojas differentiation, FIGURE 4 | Microbial community composition of the Pozas Rojas system. Relative abundance of 16S reads annotated at the genus level. Only the top 20 most abundant genera are shown in legend. Frontiers in Microbiology | www.frontiersin.org 6 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities FIGURE 5 | Non-metric multidimensional scaling plots of relative abundances at different taxonomic levels. Taxonomic categories are separated by panels (A): phylum, (B): class, (C): order, (D): family, and (E): genus. Samples are colored as follows: Pozas Rojas, purple; Lake St. Claire, pink; cold rivers of the world, red; deep water from Iowa, blue; Mediterranean Sea, yellow; Churince, black. potentiated by founder effects on each poza and their further be at least an order of magnitude lower, as is the case isolation by distance. of Svalbard subglacial ice with 700 OTUs in an equivalent Although in this study we are exploring only the water sampling coverage of three glaciers (Perini et al., 2019), community, the overall diversity in Pozas Rojas is as high as or 811 OTUs observed in five sites along the Ganjiang Churince (>5,000 OTUs, Souza et al., 2018). This is remarkable, River in China (Wang et al., 2016). It is most interesting given the small scale of the study (circa 700 × 500 m) and that despite the small scale of our study, the richness the subsampling due to the inherent sequencing technical observed here is double the species obtained in the TARA deficiencies. Other aquatic diversities have been reported to ocean survey for superficial waters and comparable to the Frontiers in Microbiology | www.frontiersin.org 7 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities FIGURE 6 | Stacked bar plot of relative abundance of phyla from global samples. Sites to which samples belong are indicated on the X-axis. observed in mesopelagic waters at the scale of the world thus supporting the marine origin of local communities (Souza (Sunagawa et al., 2015). et al., 2006). However, the observed pronounced differentiation It is possible that the large diversity found in Pozas Rojas is of the pozas from LH may be overestimated, given the acute in part the product of a more balanced stoichiometry, reducing differences in-depth between the two types of water bodies. It the environmental filtering of low phosphorus content typical is expected to get a more complete community sampling from a from CCB (Vázquez-Rosas-Landa et al., 2020). Another possible non-stratified shallow poza than a 10-m-deep lagoon from which factor is the pozas differentiation, where a particular founder the sample was taken from the surface. Hence, sampling bias effect and a cohesive community, the Eltonian niche, isolate each led by the stratification of LH most likely played a role in the poza from the next one, generating an equitable distribution, as observed divergence. well as a large diversity in most of the pozas. Another role in We observe among this large diversity a core of 150 genera such diversity and differentiation should be assigned to viruses (Figure 3), even if there is not a core at the OTU level, that we know are particularly diverse (Taboada et al., 2018) as suggesting microdiversification processes separating each pond. their constant pressure tends to increase diversity as a escape The shallow ponds are characterized by different proportions mechanism to predators. of genus from the Microbacteriaceae family with strains related It seems that each poza had a unique process of establishment to Candidatus Rhodoluna, which has been described to have [a.k.a. founder effect (Peter and Slatkin, 2015)]. In fact, after the photorhodopsin that allows for autotrophy; from the same homogenization due to the hurricane’s high waters, only pozas family are the marine genus Candidatus Aquiluna, the soil S02 and S05 share similar communities (Figure 2), even though bacteria Agrococcus, and the aquatic Leifsonia that has been they are not so close, nor part of the same environmental cluster described in continental waters (Reddy et al., 2003; Zhang (Figure 1). If we take away the lagoon (LH1), most of the unique et al., 2010; Chuon et al., 2021). Beyond the core, at the genus OTUs go away (1334) (Figure 3). Such a large number of unique level, 25 genera explain the differences between pozas with a OTUs belong mostly to 33 unique genera, with Synechococcus p > 0.015 (Supplementary Table 3). It is interesting that most being particularly abundant in LH, which is one of the most of these bacteria are heterotrophs, but only few of them are abundant Cyanobacteria in the ocean along with Prochlorococcus primary producers, as the marine Cyanobacteria Acaryochloris, (Sunagawa et al., 2015). This result is interesting but not new, an halophile that has low light chlorophyll; Aphanothece and as multiple marine lineages such as Exiguobacterium (Rebollar Prochlorothrix, common in freshwater systems; Leptolyngbya, et al., 2012), Bacillus (Alcaraz et al., 2008; Moreno-Letelier et al., found in the Arctic Ocean; and Chlorobaculum, a green sulfur 2012), and marine phages (Desnues et al., 2008) have been bacteria with an anoxygenic phototrophic metabolism. The previously found on the water systems from Cuatro Cienegas, chemoautotroph Acidithiomicrobium, an Fe-oxidizing acidophile Frontiers in Microbiology | www.frontiersin.org 8 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities and thermophile, was also found to be important in the AUTHOR CONTRIBUTIONS differentiation of the pozas, and as only few species have been described (Norris et al., 2011), it is possible that this one in MG-U: initial drafting, data analysis, figure generation, Pozas Rojas is a new species. Some genera are typical from deep discussion, literature review, and draft correction. VS and LEE: aquifers or marine environments, whereas others are typical from project coordination, discussion, and draft correction. DE-H: freshwater and soil, therefore representing potential migrants raw data processing, discussion, and literature review. JS-P and from the local environments that got “dragged” along the flood MV: initial drafting, data analysis, figure generation, discussion, water to the Pozas Rojas system. Interestingly, such migration did and literature review. LE-A: project coordination, data analysis, not seem to occur to the larger lagoon. discussion, and draft correction. MM-R: data analysis, figure When comparing the Pozas Rojas system with other water generation, discussion, and literature review. MN-M, JR-P, and systems of the world that have comparable data in MG-RAST, we CM-G: data analysis, discussion, and literature review. DM-T: observe that each site has a “signature” even at the phylum level data analysis, discussion, literature review, and draft correction. determined mostly by their most abundant groups (Figure 6). In MR-B: DNA extraction and raw data processing. MV-R-L: this case, Churince and Pozas Rojas appear similar to the cold raw data processing, data analysis, and discussion. All authors rivers of the world and phreatic mantle from Iowa, meaning, contributed to the article and approved the submitted version. continental waters. On the other hand, LH from Pozas Rojas is more similar to the epipelagic samples from the Mediterranean Sea, in particular, on their high abundance of Cyanobacteria FUNDING (Figure 6). This similarity is maintained in the NMDS plot at the phylum, order, and genus taxonomic levels only with This work was started with Fundación Carlos Slim funds and the epipelagic fraction of the Mediterranean Sea, which is continued with PAPIIT grant IG200319. This manuscript was the zone with the highest abundance of Cyanobacteria in the written with Chile ANID-CEQUA R20F0009 to VS and LEE as ocean (Figure 5). part of a larger comparison with ocean waters in the Southern tip Even though more in-depth studies are needed in this diverse of the American continent. and heterogeneous site, it is captivating that the deep aquifer possibly connected to the deep lagoon is related to the ocean, and its community structure has persisted for at least 30 million ACKNOWLEDGMENTS years that the highland of central Mexico along with the Sierra Madre Oriental uplifted isolating CCB from the Western Seaway, We would like to thank APFF of Cuatro Cienegas for their confirming that CCB is indeed a dynamic lost world where valuable help during fieldwork and sampling collections, Jerome ancestral residents mingle with more recent migrants without Verleyen and Veronica Jimenez-Jacinto for the technical support losing its particular diverse marine signature. In the last years, during bioinformatics analysis. We deeply thank Erika Aguirre efforts toward saving this amazing wetland have yielded fruit Planter and Rosalinda Tapia for support with technical logistics as 40 Ha of the wetland is in recovery as two channels that along the study. We appreciate the assistance of Felipe García diverted water for irrigation have been partially closed, giving us and Rodrigo Vázquez (both from CIECO, now Instituto the luxury of a little more time to understand the dynamics of this de Investigaciones en Ecosistemas y Sustentabilidad-UNAM) unique deep aquifer. on nutrient determinations. We also thank the Instituto de Biotecnología-UNAM for giving us access to its computer cluster and Instituto de Geofísica-UNAM for mineral determination. DATA AVAILABILITY STATEMENT SUPPLEMENTARY MATERIAL The datasets presented in this study can be found in online repositories. The names of the repository/repositories and The Supplementary Material for this article can be found accession numbers can be found below: NCBI:PRJNA785576, online at: https://www.frontiersin.org/articles/10.3389/fmicb. MG-RAST:mgp94066. 2022.825167/full#supplementary-material REFERENCES classification of small and large subunit rRNA in metagenomic data. Mol. Ecol. Res. 15, 1403–1414. doi: 10.1111/1755-0998.12399 Alcaraz, L. D., Olmedo, G., Bonilla, G., Cerritos, R., Hernández, G., Cruz, A., et al. Bonilla-Rosso, G., Peimbert, M., Alcaraz, L. D., Hernández, I., Eguiarte, L. E., (2008). The genome of Bacillus coahuilensis reveals adaptations essential for Olmedo-Alvarez, G., et al. (2012). Comparative metagenomics of two microbial survival in the relic of an ancient marine environment. Proc. Natl. Acad. Sci. mats at Cuatro Ciénegas Basin II: community structure and composition in 105, 5803–5808. doi: 10.1073/pnas.0800981105 oligotrophic environments. Astrobiology 12, 659–673. doi: 10.1089/ast.2011. Avitia, M., Escalante, A. E., Rebollar, E. A., Moreno-Letelier, A., Eguiarte, L. E., 0724 and Souza, V. (2014). Population expansions shared among coexisting bacterial Bradstreet, R. B. (1954). Kjeldahl method for organic nitrogen. Anal. Chem. 26, lineages are revealed by genetic evidence. PeerJ 2:e696. doi: 10.7717/peerj.696 185–187. Bengtsson-Palme, J., Hartmann, M., Eriksson, K. M., Pal, C., Thorell, K., Larsson, Chuon, K., Kim, S. Y., Meas, S., Shim, J-G, Cho, S.-G., Kang, K.-W., et al. D. G. J., et al. (2015). METAXA2: improved identification and taxonomic (2021). Assembly of Natively Synthesized Dual Chromophores Into Functional Frontiers in Microbiology | www.frontiersin.org 9 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities Actinorhodopsin. Front. Microbiol. 12:652328. doi: 10.3389/fmicb.2021.65 response to long-term environmental stress. Antonie Van Leeuwenhoek 104, 2328 1159–1175. doi: 10.1007/s10482-013-0038-7 Desnues, C., Rodriguez-Brito, B., Rayhawk, S., Kelley, S., Tran, T., Haynes, M., et al. Pajares, S., Souza, V., and Eguiarte, L. E. (2015). Multivariate and phylogenetic (2008). Biodiversity and biogeography of phages in modern stromatolites and analyses assessing the response of bacterial mat communities from an thrombolites. Nature 452, 340–343. doi: 10.1038/nature06735 ancient oligotrophic aquatic ecosystem to different scenarios of long-term Dick, W. A., and Tabatabai, M. A. (1977). Determination of orthophosphate environmental disturbance. PLoS One 10:e0119741. doi: 10.1371/journal.pone. in aqueous solutions containing labile organic and inorganic 0119741 phosphorus compounds. J. Environ. Q. 6, 82–85. doi: 10.2134/jeq1977. Pascual-García, A., and Bell, T. (2020). Community-level signatures of ecological 00472425000600010018x succession in natural bacterial communities. Nat. Commun. 11:2386. doi: 10. Dixon, P. (2003). VEGAN, a package of R functions for community ecology. J. Veg. 1038/s41467-020-16011-3 Sci. 14, 927–930. doi: 10.1111/j.1654-1103.2003.tb02228.x Paulson, J. N., Pop, M., and Bravo, H. C. (2013). MetagenomeSeq: statistical Elton, C. (1927). Animal Ecology. Sidgwick and Jackson. (London: The University analysis for sparse high-throughput sequencing. Biocon. Pack. 1:191. of Chicago Press). Peimbert, M., Alcaraz, L. D., Bonilla-Rosso, G., Olmedo-Alvarez, G., García-Oliva, Escobar-Zepeda, A., Godoy-Lozano, E. E., Raggi, L., Segovia, L., Merino, E., F., Segovia, L., et al. (2012). Comparative metagenomics of two microbial mats Gutiérrez-Rios, R. M., et al. (2018). Analysis of sequencing strategies and tools at Cuatro Ciénegas Basin I: ancient lessons on how to cope with an environment for taxonomic annotation: defining standards for progressive metagenomics. under severe nutrient stress. Astrobiology 12, 648–658. doi: 10.1089/ast.2011. Sci. Rep. 8:12034. 0694 García-Ulloa, M. II, Escalante, A. E., Moreno-Letelier, A., Eguiarte, L. E., and Souza, Perini, L., Gostinčar, C., and Gunde-Cimerman, N. (2019). Fungal and bacterial V. (2021). Evolutionary Rescue of an Environmental Pseudomonas otitidis in diversity of Svalbard subglacial ice. Sci. Rep. 9:20230. doi: 10.1038/s41598-019- Response to Anthropogenic Perturbation. Front. Microbiol. 11:563885. doi: 10. 56290-5 3389/fmicb.2020.563885 Peter, B. M., and Slatkin, M. (2015). The effective founder effect in a spatially García-Ulloa, M., Ponce-Soto, G. Y., González-Valdez, A., González-Pedrajo, B., expanding population. Evolution 69, 721–734. doi: 10.1111/evo.12609 Díaz-Guerrero, M., Souza, V., et al. (2019). Two Pseudomonas aeruginosa clonal Rebollar, E. A., Avitia, M., Eguiarte, L. E., González-González, A., Mora, L., Bonilla- groups belonging to the PA14 clade are indigenous to the Churince system Rosso, G., et al. (2012). Water–sediment niche differentiation in ancient marine in Cuatro Ciénegas Coahuila, México. Environ. Microbiol. 21, 2964–2976. doi: lineages of Exiguobacterium endemic to the Cuatro Cienegas Basin. Environ. 10.1111/1462-2920.14692 Microbiol. 14, 2323–2333. doi: 10.1111/j.1462-2920.2012.02784.x Grinnell, J. (1904). The origin and distribution of the chest-nut-backed chickadee. Reddy, G. S. N., Prakash, J. S. S., Srinivas, R., Matsumoto, G. I., and Shivaji, S. Auk 21, 364–382. doi: 10.2307/4070199 (2003). Leifsonia rubra sp. nov. and Leifsonia aurea sp. nov., psychrophiles from Hollister, E. B., Engledow, A. S., Hammett, A. J. M., Provin, T. L., Wilkinson, a pond in Antarctica. J. Med. Microbiol. 53, 977–984. doi: 10.1099/ijs.0.02396-0 H. H., and Gentry, T. J. (2010). Shifts in microbial community structure along Ren, L., Jeppesen, E., He, D., Wang, J., Liboriussen, L., Xing, P., et al. (2015). pH an ecological gradient of hypersaline soils and sediments. ISME J. 4, 829–838. influences the importance of niche-related and neutral processes in lacustrine doi: 10.1038/ismej.2010.3 bacterioplankton assembly. Appl. Environ. Microbiol. 81, 3104–3114. doi: 10. Junker, R. R., Lechleitner, M. H., Kuppler, J., and Ohler, L. M. (2019). 1128/AEM.04042-14 Interconnectedness of the Grinnellian and Eltonian niche in regional and local Rosado, B. H., Figueiredo, M. S., de Mattos, E. A., and Grelle, C. E. (2016). Eltonian plant-pollinator communities. Front. Plant Sci. 10:1371. doi: 10.3389/fpls.2019. shortfall due to the Grinnellian view: functional ecology between the mismatch 01371 of niche concepts. Ecography 39, 1034–1041. doi: 10.1111/ecog.01678 Kaufman, L., and Rousseeuw, P. J. (2009). Finding Groups in Data: An Introduction Rudi, K., Skulberg, O. M., Larsen, F., and Jakobsen, K. S. (1997). Strain to Cluster Analysis. (New York, NY: John Wiley & Sons), 344. characterization and classification of oxyphotobacteria in clone cultures on the Keegan, K. P., Glass, E. M., and Meyer, F. (2016). “MG-RAST, a metagenomics basis of 16S rRNA sequences from the variable regions V6, V7, and V8. Appl. service for analysis of microbial community structure and function,” in Environ. Microbiol. 63, 2593–2599. doi: 10.1128/aem.63.7.2593-2599.1997 Microbial Environmental Genomics (MEG), ed. N. J. Clifton (New York, NY: Soberón, J. (2007). Grinnellian and Eltonian niches and geographic distributions Humana Press), 207–233. doi: 10.1007/978-1-4939-3369-3_13 of species. Ecol. Lett. 10, 1115–1123. doi: 10.1111/j.1461-0248.2007.01107.x Larson, E. R., Olden, J. D., and Usio, N. (2010). Decoupled conservatism of Souza, V., Espinosa-Asuar, L., Escalante, A. E., Eguiarte, L. E., Farmer, J., Forney, Grinnellian and Eltonian niches in an invasive arthropod. Ecosphere 1, 1–13. L., et al. (2006). An endangered oasis of aquatic microbial biodiversity in the doi: 10.1890/es10-00053.1 Chihuahuan desert. Proc. Natl. Acad. Sci. 103, 6565–6570. doi: 10.1073/pnas. Lex, A., Gehlenborg, N., Strobelt, H., Vuillemot, R., and Pfister, H. (2014). UpSet: 0601434103 visualization of intersecting sets. IEEE Transact. Visual. Comput. Grap. 20, Souza, V., Moreno-Letelier, A., Travisano, M., Alcaraz, L. D., Olmedo, G., and 1983–1992. doi: 10.1109/TVCG.2014.2346248 Eguiarte, L. E. (2018). The lost world of Cuatro Cienegas Basin, a relictual Magoč, T., and Salzberg, S. L. (2011). FLASH: fast length adjustment of short reads bacterial niche in a desert oasis. Elife 7:e38278. doi: 10.7554/eLife.38278 to improve genome assemblies. Bioinformatics 27, 2957–2963. doi: 10.1093/ Su, X., Pan, W., Song, B., Xu, J., and Ning, K. (2014). Parallel-META 2.0: enhanced bioinformatics/btr507 metagenomic data analysis with functional annotation, high performance McMurdie, P. J., and Holmes, S. (2013). phyloseq: an R package for reproducible computing and advanced visualization. PLoS One 9:e89323. doi: 10.1371/ interactive analysis and graphics of microbiome census data. PLoS One journal.pone.0089323 8:e61217. doi: 10.1371/journal.pone.0061217 Sunagawa, S., Coelho, L. P., Chaffron, S., Kultima, J. R., Labadie, K., Salazar, G., Moreno-Letelier, A., Olmedo-Alvarez, G., Eguiarte, L. E., and Souza, V. (2012). et al. (2015). Structure and function of the global ocean microbiome. Science Divergence and phylogeny of Firmicutes from the Cuatro Ciénegas Basin, 348:1261359 doi: 10.1126/science.1261359 Mexico: a window to an ancient ocean. Astrobiology 12, 674–684. doi: 10.1089/ Taboada, B., Isa, P., Gutiérrez-Escolano, A. L., del Ángel, R. M., Ludert, J. E., ast.2011.0685 et al. (2018). The viral structure of the Cuatro Ciénegas Basin, a unique oasis Norris, P. R., Davis-Belmar, C. S., Brown, C. F., and Calvo-Bado, L. A. in Northern Mexico, reveals a highly diverse population at a small geographic (2011). Autotrophic, sulfur-oxidizing actinobacteria in acidic environments. scale. Appl. Environ. Microbiol. 84, e00465–18. doi: 10.1128/AEM.00465-18 Extremophiles 15, 155–163. doi: 10.1007/s00792-011-0358-3 Turner, S., Pryer, K. M., Miao, V. P., and Palmer, J. D. (1999). Investigating deep Ortiz-Álvarez, R., Cáliz, J., Camarero, L., and Casamayor, E. O. (2020). phylogenetic relationships among cyanobacteria and plastids by small subunit Regional community assembly drivers and microbial environmental rRNA sequence analysis 1. J. Eukaryot. Microbiol. 46, 327–338. doi: 10.1111/j. sources shaping bacterioplankton in an alpine lacustrine district (Pyrenees, 1550-7408.1999.tb04612.x Spain). Environ. Microbiol. 22, 297–309. doi: 10.1111/1462-2920.14 Uritskiy, G., Getsin, S., Munn, A., Gomez-Silva, B., Davila, A., Glass, B., et al. 848 (2019). Halophilic microbial community compositional shift after a rare rainfall Pajares, S., Eguiarte, L. E., Bonilla-Rosso, G., and Souza, V. (2013). Drastic changes in the Atacama Desert. ISME J. 13, 2737–2749. doi: 10.1038/s41396-019-0 in aquatic bacterial populations from the Cuatro Cienegas Basin (Mexico) in 468-y Frontiers in Microbiology | www.frontiersin.org 10 April 2022 | Volume 13 | Article 825167 García-Ulloa et al. Differentiation of Aquatic Bacterial Communities Vázquez-Rosas-Landa, M., Ponce-Soto, G. Y., Aguirre-Liguori, J. A., Thakur, Zhang, J. Y., Liu, X. Y., and Liu, S. J. (2010). Agrococcus terreus sp. S., Scheinvar, E., Barrera-Redondo, J., et al. (2020). Population genomics of nov. and Micrococcus terreus sp. nov., isolated from forest soil. Vibrionaceae isolated from an endangered oasis reveals local adaptation after Int. J. Syst. Evol. Microbiol. 60, 1897–1903. doi: 10.1099/ijs.0.013 an environmental perturbation. BMC Genom. 21:418. doi: 10.1186/s12864-020- 235-0 06829-y Vázquez-Rosas-Landa, M., Ponce-Soto, G. Y., Eguiarte, L. E., and Souza, V. (2017). Conflict of Interest: The authors declare that the research was conducted in the Comparative genomics of free-living Gammaproteobacteria: pathogenesis- absence of any commercial or financial relationships that could be construed as a related genes or interaction-related genes? Pathog. Dis. 75:ftx059. doi: 10.1093/ potential conflict of interest. femspd/ftx059 Wang, P., Chen, B., Yuan, R., Li, C., and Li, Y. (2016). Characteristics of aquatic Publisher’s Note: All claims expressed in this article are solely those of the authors bacterial community and the influencing factors in an urban river. Sci. Total and do not necessarily represent those of their affiliated organizations, or those of Environ. 569, 382–389. doi: 10.1016/j.scitotenv.2016.06.130 the publisher, the editors and the reviewers. Any product that may be evaluated in Wei, T., and Simko, V. (2017). R package “corrplot”: Visualization of a Correlation this article, or claim that may be made by its manufacturer, is not guaranteed or Matrix (Version 0.84). Available online at: https://github.com/taiyun/corrplot endorsed by the publisher. (accessed August 2021). Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. (New York, NY: Copyright © 2022 García-Ulloa, Souza, Esquivel-Hernández, Sánchez-Pérez, Springer-Verlag). Espinosa-Asuar, Viladomat, Marroquín-Rodríguez, Navarro-Miranda, Ruiz- Wolaver, B. D., Crossey, L. J., Karlstrom, K. E., Banner, J. L., Cardenas, M. B., Ojeda, Padilla, Monroy-Guzmán, Madrigal-Trejo, Rosas-Barrera, Vázquez-Rosas-Landa C. G., et al. (2013). Identifying origins of and pathways for spring waters in and Eguiarte. This is an open-access article distributed under the terms of a semiarid basin using He. Sr, and C isotopes: Cuatrociénegas Basin, Mexico. the Creative Commons Attribution License (CC BY). The use, distribution Geosphere 9, 113–125. doi: 10.1130/ges00849.1 or reproduction in other forums is permitted, provided the original author(s) Yan, Q., Stegen, J. C., Yu, Y., Deng, Y., Li, X., Wu, S., et al. (2017). Nearly a decade- and the copyright owner(s) are credited and that the original publication long repeatable seasonal diversity patterns of bacterioplankton communities in this journal is cited, in accordance with accepted academic practice. No in the eutrophic Lake Donghu (Wuhan, China). Mol. Ecol. 26, 3839–3850. use, distribution or reproduction is permitted which does not comply with doi: 10.1111/mec.14151 these terms. Frontiers in Microbiology | www.frontiersin.org 11 April 2022 | Volume 13 | Article 825167