Metabarcoding of the phytotelmata of Pseudalcantarea grandis (Bromeliaceae) from an arid zone

Background Pseudalcantarea grandis (Schltdl.) Pinzón & Barfuss is a tank bromeliad that grows on cliffs in the southernmost portion of the Chihuahuan desert. Phytotelmata are water bodies formed by plants that function as micro-ecosystems where bacteria, algae, protists, insects, fungi, and some vertebrates can develop. We hypothesized that the bacterial diversity contained in the phytotelma formed in a bromeliad from an arid zone would differ in sites with and without surrounding vegetation. Our study aimed to characterize the bacterial composition and putative metabolic functions in P. grandis phytotelmata collected in vegetated and non-vegetated sites. Methods Water from 10 individuals was sampled. Five individuals had abundant surrounding vegetation, and five had little or no vegetation. We extracted DNA and amplified seven hypervariable regions of the 16S gene (V2, V4, V8, V3–6, 7–9). Metabarcoding sequencing was performed on the Ion Torrent PGM platform. Taxonomic identity was assigned by the binning reads and coverage between hit and query from the reference database of at least 90%. Putative metabolic functions of the bacterial families were assigned mainly using the FAPROTAX database. The dominance patterns in each site were visualized with rank/abundance curves using the number of Operational Taxonomic Units (OTUs) per family. A percentage similarity analysis (SIMPER) was used to estimate dissimilarity between the sites. Relationships among bacterial families (identified by the dominance analysis and SIMPER), sites, and their respective putative functions were analyzed with shade plots. Results A total of 1.5 million useful bacterial sequences were obtained. Sequences were clustered into OTUs, and taxonomic assignment was conducted using BLAST in the Greengenes databases. Bacterial diversity was 23 phyla, 52 classes, 98 orders, 218 families, and 297 genera. Proteobacteria (37%), Actinobacteria (19%), and Firmicutes (15%) comprised the highest percentage (71%). There was a 68.3% similarity between the two sites at family level, with 149 families shared. Aerobic chemoheterotrophy and fermentation were the main metabolic functions in both sites, followed by ureolysis, nitrate reduction, aromatic compound degradation, and nitrogen fixation. The dominant bacteria shared most of the metabolic functions between sites. Some functions were recorded for one site only and were related to families with the lowest OTUs richness. Bacterial diversity in the P. grandis tanks included dominant phyla and families present at low percentage that could be considered part of a rare biosphere. A rare biosphere can form genetic reservoirs, the local abundance of which depends on external abiotic and biotic factors, while their interactions could favor micro-ecosystem resilience and resistance.


INTRODUCTION
Phytotelmata are water bodies formed by plants that function as micro-ecosystems (Benzing, 2000). The community comprises bacteria, cyanobacteria, protists, fungi, green algae, mosses, vascular plants, insects, crustaceans, and a few vertebrates (Benzing, 2000;Kitching, 2001;Brandt, Martinson & Conrad, 2016). Under natural conditions, organisms are frequently replaced, and the system has been used as a study model for food webs (Mogi, 2004). Phytotelmata are most frequently found in tropical areas but can also occur in temperate forests, swamps, and deserts. In arid environments, the phytotelmaassociated micro-ecosystem is defined by the seasonality of water availability. Once water accumulates following the rains, growth occurs in the aquatic biota that is well adapted to temporary environments, significantly increasing the diversity of aquatic organisms in the area (Calhoun et al., 2017).
Although different plant families form phytotelmata, the Bromeliaceae have various anatomical, morphological, and physiological adaptations that allow them to grow in areas with wide resource variations (Giongo et al., 2019). For example, the leaves are arranged in a tight rosette, and the plant epidermis is covered with trichomes that absorb humidity and nutrients, allowing the plants to grow in arid environments with scarce nutrients (Benzing, 2000;Goffredi, Kantor & Woodside, 2011). Pseudalcantarea grandis (Schltdl.) Pinzón & Barfuss is a bromeliad found in saxicolous habitats, up to 2.5 m in height and with a branched inflorescence present in March and April. It is native to centraleastern Mexico to Honduras (Rzedowski, 2006). The species thrives on canyon cliffs of the major rivers of the northeastern Bajío region, Mexico, at altitudes ranging from 400 to 1,600 m asl. Due to the inaccessibility of its populations, it presents no particular conservation problems.
Characterization and identification of organisms contained in environmental samples can be conducted using different approaches, such as sample culture, target sequencing, metabarcoding, metatranscriptomics, and metagenomics. The diversity of specific groups in the tank bromeliads has been analyzed with targeted sequencing on ciliates and vertebrates (Brozio et al., 2017;Simão et al., 2017). Using metatranscriptomics, Goffredi, Jang & Haroon (2015) found 450 species of Archaea and bacteria in Vriesea platynema Gaudich. (Bromeliaceae) tanks. Metabarcoding is the direct analysis of DNA fragments contained in an environmental sample (Cabral et al., 2018). This technique allows the identification of microorganisms with no need for culturing (Rodríguez-Nuñez, Rullan-Cardec & Rios-Velazquez, 2018). Metabarcoding has been used to identify bacterial and eukaryotic biodiversity in the phytotelmata of Sarracenia purpurea L. (Sarrraceniaceae) (Grothjan & Young, 2019). In tank bromeliads, bacterial metabarcoding has been used in five studies, four in Brazil and one in Puerto Rico, all in tropical forests Louca et al., 2017;Simão et al., 2017;Rodríguez-Nuñez, Rullan-Cardec & Rios-Velazquez, 2018;Giongo et al., 2019;Simão et al., 2020). To our knowledge, however, arid zone bromeliads have not been studied.
The biotic composition of the phytotelmata depends on the species, its location, and local factors that affect water conditions (Benzing, 2000;Louca et al., 2017;Males, 2016). Bromeliad tanks form a unique freshwater environment that differs in oxygen concentration and pH from the external environment, thus providing a habitat for a diverse community (Goffredi, Kantor & Woodside, 2011). The phytotelmata in bromeliads from tropical forests can contain methanogens, which are microorganisms responsible for carbon cycling (Goffredi, Kantor & Woodside, 2011). When comparing the community of archaea and methanogens in phytotelmata from different tank water volume, it was found that the methane cycle formation in the phytotelma decreases during dry periods in neotropical forests (Brandt, Martinson & Conrad, 2016). Identifying the bacterial communities of bromeliad phytotelmata from different ecological niches can help to understand their interaction with the metabolism of the host plant (Louca et al., 2017). The phytotelmata of P. grandis constitutes a temporary aquatic ecosystem in a desert, and its biodiversity has not been studied. Although water availability is highly seasonal, we hypothesized that the tank bacterial composition will differ in sites with and without surrounding vegetation. Our study aimed to characterize the bacterial composition and putative metabolic functions in P. grandis phytotelmata collected in vegetated and non-vegetated sites.

Site descriptions, plant selection, and sampling
The study site is located in the Las Angosturas canyon, also known as Barranca Tolimán,in Zimapán,Hidalgo,in central Mexico (20 50.933′N,99 26.7′W, 900 masl) (Fig. 1). The area is located in the southernmost portion of the Chihuahuan desert (Hernández & Gómez-Hinostrosa, 2005) and constitutes a local floristic region of high endemism (Medellín-Leal, 1982). The exact location of the study area does not feature in any geomorphological or geological publications. However, adjacent canyons in the same region have been subjected to detailed studies (Segestrom, 1961;Carrillo, 1981;Carrillo & Sutter, 1981;Arévalo, 1991). The geological formations are Trancas (Late Jurassic, Early Cretaceous), el Doctor (Middle Cretaceous), and Soyatal (Upper Cretaceous), formed by a combination of calcareous rocks alternated with calcareous limestones, calcareous lutites, and sandstones. Structurally, the canyon is formed by rocky vertical cliffs at 80-90 angles. The P. grandis plants grow on sandstone rocks on the vertical cliffs ( Fig. 2) of the El Doctor formation. Ten individuals of 50 cm or more in diameter were sampled on cliffs: five with little or no vegetation ( Fig. 2A) and five with abundant surrounding vegetation (Fig. 2B). Four of our vegetated sample sites had a NE orientation and one a NW orientation; all sites were surrounded by either xerophytic scrub or tropical deciduous forest. The non-vegetated sites all faced N. The plant species surrounding the sample sites were identified and recorded (Table 1). Water samples were collected in June 2018 during the rainy season since the plants are dry for the rest of the year, either empty or full of debris (Figs. 2C, 2D). Experiments were approved by the "Comité de Bioética de la Facultad de Ciencias Naturales" bioethics committee (39FCN2019). Bromeliads were reached by rappel (Figs. 2E, 2F). NestÒ cell scrapers were used to scratch the inside of each tank, and the water in the bromeliad was vigorously shaken in order to obtain a homogeneous sample. Water volumes of 50 to 100 ml were collected using 10 ml sterile serological pipettes. Samples were stored in 50 ml conical Falcon tubes, transported on dry ice, and stored at −79 C until processed.

DNA extraction and sequencing
The five samples of each site were homogenized and pooled. A total of 100 ml of sampled water was filtered through a 0.22 µm nitrocellulose MilliporeÒ membrane. The membrane was then frozen and macerated in liquid nitrogen. DNA was extracted in triplicate with the QIAmp DNA extractionÒ kit following the manufacturer's instructions. DNA quality and quantity were evaluated using spectrophotometry in a NanoDropÒ instrument. PCR amplicons of seven hypervariable regions of the 16S gene were amplified with two primer sets, the first targeting V2, V4, V8, and the second V3-6, 7-9, with the Ion 16S TM Metagenomics kit (Thermo Fisher Scientific, Waltham, MA, USA), following the manufacturer's protocol. The metabarcoding sequencing was performed on the Ion Torrent PGM platform, and the amplicons were purified with AgencourtÒ AMPureÒ XP. The Ion Plus Fragment Library kit protocol was followed in order to construct the libraries. Fragment presence, size, and concentration were analyzed using a Bioanalyzer 2100 with the High Sensitivity DNA assay (Agilent, Santa Clara, CA, USA). Libraries were quantified using real-time PCR to obtain an equimolar dilution factor for mixing the libraries. Templates were prepared via an emulsion PCR in the Ion One Touch System (Life Technologies, Carlsbad, CA, USA) and quantified in a fluorometer in QubitÒ 3.0 (Thermo Fisher Scientific, Waltham, MA, USA). The template was loaded in the PGM 318 TM chip using the sequencing kit for 400 base pairs, following the Ion PGM TM Hi-Q TM View Sequencing Kit protocol. Data analysis

Bioinformatic analysis
Bacteria were determined using Ion Reporter TM . Sequencing results were analyzed using the metagenomics application for multiple groups based on the Greengenes v13.5 database. Primers used for amplification were identified, and a minimum sequence length of 150 bp was defined. To assign taxonomic identity, we considered two criteria: the binning reads had to be repeated at least 10 times, and the coverage between hit and query from the reference database had to be at least 90%.

Analysis of bacterial composition between sites
Bacterial families were ordered by taxonomic hierarchy for each site, and a richness stacked barplot was produced at order and family level with Microsoft Excel tools. The bacterial composition of the two sites vegetated (V) and non-vegetated (NV) was compared using the Sørensen similarity coefficient based on a presence/absence matrix for bacterial families, and a Venn diagram was generated using vegan and VennDiagram packages in R Studio v3.6.1 (R Core Team, 2019). In addition, the dominance patterns of bacterial families were visualized with rank/abundance curves, using the number of OTUs per family. A percentage similarity analysis (SIMPER) was used to estimate the dissimilarity between sites. SIMPER was performed with the composition and number of OTUs per family, a data pretreatment by square root-transformation, and the Bray-Curtis similarity coefficient. A shade plot was constructed using the most important bacterial families, according to dominance and contribution to the dissimilarity between the vegetated and non-vegetated sites. These bacterial families were selected with the bacterial dominance analysis and SIMPER results, considering a cumulative contribution of~40% in both. In this shade plot, a matrix of family composition and number of OTUs was used. For the classification of families, a Whittaker association coefficient was used with data previously standardized to percentages, and the group average linkage method. In the samples from vegetated and non-vegetated sites, a Bray-Curtis similarity and a square-root transformation were used.

Metabolic functions
To identify putative metabolic functions, we used the FAPROTAX database v.1.2.4 . This database assigns a putative metabolic function to each OTU based on the literature and, for some taxa, associates this function with cultured taxa with a verified function in the same taxonomic group. The current bacterial diversity not recognized under culture is high, and therefore the generalized assignment may change in future studies. However, this database provides information on 4,600 taxa (Louca et al., 2017). We analyzed the data in two ways: combined and separated (vegetated and non-vegetated sites), and with data from each site separately. Function was assigned at the family and genus level whenever possible. The putative taxa function that was absent from the FAPROTAX database was inferred based on the available literature. We looked for the family name and then reviewed its metabolic functions (Bergey & Holt, 2005;Louca et al., 2017). Furthermore, the relationship between bacterial families identified by the dominance and SIMPER analysis, and their respective putative functions, was analyzed with another shade plot. This analysis was performed with a binary matrix based on the Sørensen similarity coefficient to associate families and functions using the group average linkage method. The range/abundance curves, SIMPER, and shade plots were generated in PRIMER 7 7.0.21 (Clarke & Gorley, 2015).

DISCUSSION
The results of this study indicate that the composition of bacterial families in the phytotelmata of P. grandis is similar between the vegetated and non-vegetated sites. Nevertheless, they present a different dominance pattern as a function of the richness of OTUs associated with these families. Bacterial richness in P. grandis is composed of 23 phyla and 218 families. This result contrasts with Aechmea bromeliifolia and A. nudicaulis, each of which contain 51 phyla (Rodríguez-Nuñez, Rullan-Cardec & Rios-Velazquez, 2018), and with the 30 phyla reported in Aechmea gamosepala, Vriesea friburgensis, and V. platynema (Simão et al., 2020). However, at the family level, we found a greater richness in P. grandis compared to Aechmea nudicaulis (81 families , A. nudicaulis and Neoregelia cruenta (56 families Louca et al., 2017), and A. gamosepala and V. platynema (103 families, Giongo et al., 2019). We considered that the families with a low richness of OTUs that contribute a low percentage (<1%) to the diversity of P. grandis could possibly be considered as a rare biosphere (Pedrós-Alió, 2012;Jousset et al., 2017). Although the record of these families could be a product of chance rather than ecological forces, the triplicate sequencing decreases such probability. In this study, another factor contributing to the detection of these OTUs was the use of seven hypervariable regions of the 16S. Some studies demonstrate that these regions vary in sensitivity and level of informativeness for different approaches (Yang, Wang & Qian, 2016;Fiannaca et al., 2018;Huttenhower et al., 2012;Soergel et al., 2012;D'Amore et al., 2016;Zheng et al., 2015;Chakravorty et al., 2007).
The dominant metabolic function within the bacterial community in tank bromeliad is the decomposition of complex organic compounds accumulated as vegetal detritus . Members of the Phylum Actinobacteria are saprophytes that decompose a wide spectrum of plant and animal remains (Zhang et al., 2017). They also occur in polluted environments of both terrestrial and aquatic ecosystems (Rosenberg, Delong & Thompson, 2014). Proteobacteria are the dominant group in soil microbial communities (Zhang et al., 2017), as well as in bromeliad phytotelmata Louca et al., 2017). Many Firmicutes can also decompose organic debris, resist high temperatures, and remain in dehydrated environments by inactivity (Parkes & Sass, 2009). Their presence in the P. grandis tanks suggests the occurrence of a nutrient recycling process, which provides resources for both the associated biota and the plant itself.
Despite some differences in bacterial taxonomic diversity in P. grandis between vegetated and non-vegetated sites, the dominant bacteria share most of the metabolic functions. The six main functions, aerobic chemoheterotrophy, fermentation, ureolysis, nitrate reduction, aromatic compound degradation, and nitrogen fixation, are prominent, since these are carried out by the families with a greater amount of OTUs in both sites. The first three functions mentioned above occur in equal percentages when the 212 families were included. Aechmea nudicaulis (L.) Griseb. (Bromeliaceae) and Sarracenia purpurea L. (Sarraceniaceae) present different bacterial composition in their phytotelmata, but with similar functions Grothjan & Young, 2019). However, when the geochemical conditions of the tanks of A. nudicaulis and Neoregelia cruenta (Graham) L.B. Sm. (Bromeliaceae) are compared, functional community structure is strongly correlated with the different ecological conditions provided by the vegetal cover and access to freshwater (Louca et al., 2017). In P. grandis when the total family richness is considered a slight decrease in cellulolysis (7% vs. 5%) and a slight increase in fermentation (21% vs. 23%) were detected in vegetated compared to non-vegetated sites. Moreover, some putative functions were recorded only in one site and related to families with the lowest richness of OTUs. These differences could be related to environmental factors that were not considered in this study. More studies are required to gather conclusive evidence in this regard.
The bacterial composition of P. grandis shows differences between sites in terms of the exclusive families, relative abundance of OTUs, and percentages of putative metabolic functions performed. Although the two sites share 68.3% of their composition, the unshared families suggest variations in the physiochemical conditions of the phytotelmata at each site. The bacterial community of vegetated site presents families which require an acidic pH and high levels of organic carbon and nitrogen compounds. For example, Alicyclobacillaceae grows in acid environments produced by carbohydrates (Stackebrandt, 2014). Deinococcaeae can live with high radiation levels (Murray, 1992), and Nitrosomonadaceae play significant roles in control of the nitrogen cycle in freshwater environments as ammonia oxidizers (Prosser, Head & Stein, 2014) (Fig. 4). In contrast, the sample from the non-vegetated site contained families with metabolic functions that are associated with autotrophic organisms, and others adapted to carbon and oxygen scarcity that utilize inorganic nitrogen and sulfur compounds deposited by rock sediments in their life cycle. Some of the families are Fusobacteriaceae that ferment carbohydrates and can live in anaerobic sediments (Olsen, 2014). Desulfobacteraceae are strictly anaerobic sulfate-reducing bacteria that grow best at moderate temperatures (Kuever, 2014). Desulfuromonadaceae are found in anoxic environments and are associated with methanogens and phototrophic green sulfur bacteria (Greene, 2014) (Figs. 4, 5). The families Frankiaceae, Rickettsiaceae, and Thioalkalispiraceae also perform methylotrophy (i.e. they can obtain energy from single-carbon compounds). The largest number of families belongs to the orders Actinomycetales and Rhizobiales, taxa that degrade plant debris and comprise genera (such as Streptomyces and Rhizobium) that present symbiotic relationships with plants. Their function in P. grandis is as degraders and symbionts, promoting plant growth and maintaining the ecosystem formed inside the bromeliad. The differences in the orders and families of bacteria unique to each site indicate that, when the phytotelma is exposed, the biota will mostly be autotrophic and will utilize the rock sediments from the slope (chemoautotrophs) and sunlight (phototrophs) for their metabolic functions.
The bacterial diversity found in the tank suggests that the organisms that inhabit these small aquatic microhabitats take advantage of water availability to develop. After the dry season, endospores in the tank, or from the environment around the tank (e.g., in the air, in the debris) proliferate quickly during the short rainy season and are specialized in the decomposition of complex organic compounds. Rare biosphere bacteria (OTUs or species with frequencies less than or equal to 1% (Pedrós-Alió, 2012)) play important ecological roles as drivers of ecosystem key functions. They are also considered genetic reservoirs, the abundance of which depends on external abiotic and biotic factors. Their interactions could favor micro-ecosystem resilience and resistance (Coveley, Elshahed & Youssef, 2015;Jousset et al., 2017). We found that a few families also present in low frequencies have putative metabolic functions recorded for one site only. They include Alcanivoracaceae, which are involved in aliphatic non-methane hydrocarbon degradation and oil bioremediation. However, the presence of these families should be treated with some caution. Future studies on tank bromeliads should address the relationship between the rare families and the maintenance of the micro-ecosystem.

CONCLUSIONS
We hypothesized that bacterial diversity in the phytotelmata from an arid zone would differ in sites with and without surrounding vegetation. Slight differences were found for Pseudalcantarea grandis in taxonomic richness, number of OTUs for the dominant and exclusive families, and the putative metabolic functions performed in each site. The non-vegetated site was richer in families and exclusive OTUs than the vegetated site. In the latter, families such as Deinococcaeae and Nitrosomonodaceae prefer an acidic pH and high levels of nutrients. The phytotelma of the non-vegetated site contain families such as Fusobacteriaceae and Desulfobacteraceae that thrive under carbon and oxygen shortage and can metabolize inorganic and sulfur compounds. The organisms that inhabit the small ephemeral aquatic microhabitats are well adapted to prolonged dry periods and development quickly in water presence. Their taxonomic variation could fulfill specialized functions in the degradation of organic matter, photo-or chemoautotrophy depending on the exposure to different conditions. Our study is the first to characterize the P. grandis microbiome and the information generated will be of utility to new studies in tank bromeliads and related groups.

Author Contributions
José Alan Herrera-García performed the experiments, analyzed the data, prepared figures and/or tables, field work, and approved the final draft. Mahinda Martinez conceived and designed the experiments, analyzed the data, authored or reviewed drafts of the paper, field work, and approved the final draft. Pilar Zamora-Tavares conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Ofelia Vargas-Ponce conceived and designed the experiments, performed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Luis Hernández-Sandoval analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, field work, and approved the final draft. Fabián Alejandro Rodríguez-Zaragoza analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.

Field Study Permissions
The following information was supplied relating to field study approvals (i.e., approving body and any reference numbers): Field sampling was approved by the Facultad de Ciencias Naturales bioethics committee (39FCN2019).

Data Availability
The following information was supplied regarding data availability: The sequences of the seven hypervariable regions of the 16S are available at NCBI PRJNA685432.