The polyps of A. coerulea and A. solida were originally obtained from the laboratory of Agustin Schiariti (Instituto Nacional de Investigación y Desarrollo Pesquero). A. coerulea were collected from the USA, China, Japan, Spain, and France and abbreviated as AC_USA, AC_CHI, AC_JAP, AC_SPA and AC_FRA, respectively. A. solida were collected from Israel and Slovenia and abbreviated as AS_ISR and AS_SLO, respectively, as detailed in Table 1 . The polyps of the different Aurelia populations were kept separately in tanks. All polyps were cultured at 15°C in ambient fresh seawater with a salinity of 30 practical salinity units and kept on a day:night lighting rhythm of 12:12 h for 6 months of incubation time. Freshly hatched Artemia salnia were fed to the Aurelia polyps as their sole food source once every two days during the incubation period. Homogenization of cultivation conditions was performed to ensure stability and monovariability (i.e., host genotypes) of the symbiotic microbial communities of polyps from different Aurelia species and populations. The polyps were incubated in sterile seawater (0.22-μm filtered) baths for 1 day to clear the digestive system before DNA extraction. Species were identified based on 16S mtRNA pairwise sequence alignment technology.

Bacterial DNA of Aurelia polyps from 7 populations was isolated using a Wizard genomic DNA purification kit (Promega, Madison, WI, USA) according to the manufacturer’s protocols. Each Aurelia population encompassed 5 replicates (10 polyps each) to yield read libraries sufficient for analysis. PCR amplification of the V4 hypervariable region of the 16S rRNA gene was performed with the primers 515F (5’-GTGCCAGCMGCCGCGGTAA-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’) (Hugerth et al., 2014). All PCRs were carried out in 30 µL reactions with 15 µL of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs, USA), 0.2 µM forward and reverse primers, and approximately 10 ng of template DNA. Thermal cycling consisted of initial denaturation at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s and elongation at 72°C for 30 s, with a final elongation at 72°C for 5 min. PCR products were detected by electrophoresis in a 2% (w/v) agarose gel. PCR amplicons of each sample with bright bands were mixed in equal-density ratios and purified with a GeneJET™ Gel Extraction Kit (Thermo Scientific, USA). Sequencing libraries were generated using an Ion Plus Fragment Library Kit 48 Rxns (Thermo Scientific) following the manufacturer’s recommendations. Then, the library concentration was assessed with a Qubit@ 2.0 Fluorometer (Thermo Scientific). The amplicon libraries were sequenced on the Ion S5™ XL platform at Novogene Bioinformatics Technology Co., Ltd. (Beijing, China).

Single-end reads were assigned to samples based on their unique barcode and truncated by cutting off the barcode and primer sequence. Quality filtering of the raw reads was performed under specific filtering conditions to obtain high-quality clean reads according to the Cutadapt (v1.9.1, http://cutadapt.readthedocs.io/en/stable/) quality control process (Kechin et al., 2017). The reads were compared with the reference database (Silva 132 database, https://www.arb-silva.de/) (Quast et al., 2013) using the UCHIME algorithm (UCHIME algorithm, http://www.drive5.com/usearch/manual/uchime_algo.html) (Edgar et al., 2011) to detect chimera sequences and then the chimera sequences were removed to obtain clean reads (Haas et al., 2011). The singleton and non-target sequences were removed from the analysis.

Sequence analyses were performed using UPARSE software (UPARSE v7.0.100 http://drive5.com/uparse/) (Edgar, 2013). First, sequences with ≥ 97% similarity were assigned to the same operational taxonomic unit (OTU), and a representative sequence for each OTU was screened for further annotation. Then, the Silva 132 database (https://www.arb-silva.de/) (Quast et al., 2013) was used to annotate each representative sequence with taxonomic information based on the Mothur algorithm, and a representative sequence for each OTU was assigned to a taxonomic level using the RDP classifier (Edgar, 2013). Furthermore, multiple sequence alignments were conducted using MUSCLE software (version 3.8.31, http://www.drive5.com/muscle/) (Edgar, 2004) to study the phylogenetic relationships among different OTUs and the divergence in the dominant species among different samples (groups).

Microbial communities normally consist of a few abundant and many rare species (Easson et al., 2020). In this study, the thresholds for rare, conditionally rare, and abundant taxa were defined based on relative sequence abundance cutoffs, with reference to recent publications (Liu et al., 2017). “Rare taxa” were defined as OTUs with a relative sequence abundance < 0.01% in all samples. “Conditionally rare taxa” were defined as OTUs that were rare (relative sequence abundance < 0.01%) in some but not all samples and were never abundant (relative sequence abundance ≥ 1%). “Abundant taxa” were defined as the OTUs that did not fall in either the rare or conditionally rare categories. “Core taxa” was defined as the OTUs present in all Aurelia polyps in this study.

The OTU abundance data were normalized corresponding to the sample with the fewest sequences for further analysis of alpha diversity and beta diversity. The Shannon and CHAO1 diversity indices were calculated based on the normalized OTU matrix with QIIME (version 1.7.0) and visualized with the ggplot2 package of R software (version 2.15.3). Nonparametric Kruskal–Wallis tests, nonparametric Mann–Whitney U tests, and t tests were performed to identify significant differences in functional and relative bacterial abundance between different samples. The beta diversity of samples was calculated based on Bray–Curtis dissimilarity at the OTU level and used to perform principal coordinate analysis (PCoA), which was visualized with the ggplot2 and vegan packages of R software (version 2.15.3). In addition, Wilcoxon match-pair tests and permutational analysis of molecular variance (PERMANOVA) were constructed to test for significant differences in microbial alpha and beta diversity between Aurelia polyp species using the vegan package of R software (version 2.15.3). Furthermore, the relative sequence abundance of samples at diverse classification levels was visualized by the ggplot2 and reshape2 packages of R software (version 2.15.3). The relative sequence abundance of functional bacterial communities was calculated based on the FAPROTAX database and visualized by the pheatmap package of R software (version 2.15.3).

After filtering the raw data, a total of 2,456,681 clean reads from 35 Aurelia polyp samples were obtained, with an average sequence length of 372 bp. A total of 1,213 OTUs were present in the Aurelia polyp samples, clustered into 36 phyla, 49 classes, 107 orders, 190 families, and 400 genera.

A total of 159 OTUs were present in all A. coerulea polyp groups, accounting for 82.06% (AC_SPA)–98.58% (AC_CHI) of the total relative sequence abundance ( Figure 1 ). In the A. coerulea polyps, 36 bacterial phyla, 48 classes, 98 orders, 179 families, and 366 genera were detected. The microbiomes of A. solida polyps contained 25 phyla, 40 classes, 84 orders, 145 families, and 257 genera. A total of 334 OTUs were shared between both groups of A. solida polyps and accounted for 96.64% and 98.97% of the total relative sequence abundance in AS_ISR and AS_SLO, respectively ( Figure 1 ). Each Aurelia polyp group contained unique OTUs, which all had low relative sequence abundance (< 1%, Figure 1 ).

At the phylum and class levels, bacterial communities associated with both Aurelia polyp species were dominated by Proteobacteria (mainly classes Alphaproteobacteria and Gammaproteobacteria), Firmicutes (mainly class Clostridia), Bacteroidetes, Acidobacteira and Actinobacteria, together comprising 97.01% (AC_SPA) – 99.71% (AC_JAP) of the total richness ( Figure 2 ). The relative abundances of the other 31 phyla in each Aurelia polyp population were <1%, together comprising 0.29% (AC_JAP) – 2.99% (AC_SPA) of the total richness ( Figure 2 ). At the family and genus levels, the polyp microbiomes were dominated by Sphingomonadaceae (mainly the genera Sphingomonas and Sphingobacterium), Rhizobiaceae (mainly the genus Phyllobacterium) and Burkholderiaceae (mainly the genus Ralstonia), comprising 53.2% (AC_SPA)–95.5% (AC_CHI) of the total bacterial abundance associated with A. coerulea and A. solida ( Figure 2 ).

At the family level, unidentified Clostridiales (Mann–Whitney U test, p = 0.039) were significantly more abundant in A. solida than in A. coerulea, while Sphingomonadaceae (t test, p = 0.038) and Rhizobiaceae (t test, p = 0.024) were more abundant in A. coerulea. At the genus level, the A. coerulea polyps had significantly higher relative sequence abundance of the genera Sphingomonas (t test, p = 0.043) and Phyllobacterium (t test, p = 0.023) than A. solida polyps.

To further reveal the inter- and intraspecific differences in the microbiomes of Aurelia polyps, the microbes of the two Aurelia polyp species were compared based on microbial alpha diversity and beta diversity, as detailed in Figure 3 .

Regarding interspecific differences in microbial alpha diversity, A. coerulea polyps had a lower average Shannon index (2.998) than A. solida polyps (3.119), as well as a lower average CHAO1 index (219 and 237 in A. coerulea and A. solida polyps, respectively, Figure 3 ). However, no significant difference between the A. coerulea and A. solida groups was detected in either the Shannon indices or CHAO1 indices (Wilcoxon match-pairs test, p = 0.5566, Figure 3 ). Regarding intraspecific differences, AC_CHI had the highest average Shannon indices of the A. coerulea groups (3.367), followed by AC_SPA, AC_FRA, AC_USA, and AC_JAP had the lowest (2.635). The highest average CHAO1 index of the A. coerulea polyps was found in AC_SPA (285), followed by AC_USA, AC_FRA and AC_JAP, and the lowest was found in AC_CHI (186). Within the A. solida groups, AS_ISR had a higher average CHAO1 index (263) than AS_SLO (212) but had a lower mean Shannon index (2.449 in AS_ISR and 2.831 in AS_SLO). Nonetheless, no significant difference was found between the Shannon and CHAO1 indices of polyps in the A. coerulea and A. solida groups (Wilcoxon match-pairs test, p > 0.05, Figure 3 ).

In terms of the interspecific comparison, PCoA based on Bray–Curtis distance demonstrated a significant difference in microbial beta diversity between the microbiomes of the A. coerulea and A. solida polyps (PERMANOVA test, p = 0.016, Figure 3 ). Furthermore, significant differences were identified between the abundant and core microbial taxa associated with A. coerulea polyps and A. solida polyps (PERMANOVA test, p = 0.044 and p = 0.03, respectively), but there was no significant difference between the rare microbial taxa associated with the two species (PERMANOVA test, p = 0.403, Table 2 ). Similarity percentage analysis (Simper) illustrated that the genera Stenotrophomonas, Ralstonia, and Sphingobacterium were the main contributors to the difference in the microbial community between the two species, which contributed 26.8%, 26.1%, and 10.6% to the difference value, respectively ( Figure 4 ). In terms of the intraspecific comparison, PCoA showed no significant intraspecific variation among the groups of A. coerulea polyps (PERMANOVA test, p = 0.60, Figure 3 ) or A. solida polyps (PERMANOVA test, p = 0.484, Figure 3 ). Furthermore, no significant difference was detected in the rare, conditionally rare, abundant taxa, or core taxa of the intraspecific groups (PERMANOVA test, p > 0.05, Table 2 ).

In total, 20 OTUs were found in all Aurelia polyp samples. These 20 OTUs were defined as the “core microbiome” of Aurelia polyps, accounting for 57.19% (AC_SPA)–96.86% (AC_CHI) of the total microbial richness ( Figure 5 ). The genera Sphingomonas, Ralstonia, and Phyllobacterium were prominent in the core microbiome of the Aurelia polyps, accounting for 79.76% (AS_ISR)–92.64% (AS_SLO) and 57.19% (AC_SPA)–96.13% (AC_CHI) of the total core microbial richness in A. coerulea polyps and A. solida polyps, respectively ( Figure 5 ). Notably, 6 OTU members of the core microbiome had significantly different abundances between A. coerulea and A. solida polyps: OTU_2 (Mann–Whitney U test, p = 0.034), OTU_1043 (Mann–Whitney U test, p = 0.026), and OTU_1244 (Mann–Whitney U test, p = 0.037) all belonging to the genus Ralstonia, were significantly more abundant in A. solida polyps than in A. coerulea polyps. However, OTU_3 (genus Phyllobacterium, t test, p = 0.026), OTU_1100 (genus Sphingomonas, t test, p = 0.026), and OTU_1287 (genus Sphingomonas, Mann–Whitney U test, p = 0.011) were significantly more abundant in A. coerulea polyps than in A. solida polyps.

Ecological functional annotation of polyp-associated microbial communities was conducted based on the FAPROTAX database. A total of 1,247 functional assignments for 8, 326 OTUs were obtained. OTUs without any functional annotation were excluded from the analysis. Functional group abundances in each sample were calculated as the cumulative relative sequence abundance of OTUs assigned to each functional group after normalizing by the cumulative abundances of OTUs associated with at least one function. Generally, chemoheterotrophy (42.7~70.0% of total OTUs) and anaerobic chemoheterotrophy (36.8~65.5% of total OTUs) were the primary functions in both A. coerulea and A. solida polyps, followed by functions related to the N cycle, such as nitrate reduction, nitrogen respiration, and denitrification ( Figure 6 ). Furthermore, microbial communities associated with A. coerulea polyps had significantly higher abundances of chemoheterotrophic bacterial taxa than those associated with A. solida polyps (Mann–Whitney U test, p = 0.004).

Bacterial communities associated with hosts are shaped by host selection and environmental conditions (Gould et al., 2018). Previous researchers have highlighted environmental factors shaping the structure of microbial communities in hosts such as coral (Zhang et al., 2014; Osman et al., 2020), sponge (Zhang et al., 2014; Easson et al., 2020), and sea anemone (Mortzfeld et al., 2016; Morelan et al., 2019). In this study, despite long-term cultivation in the same environment and feed, we found interspecific variation in beta diversity among Aurelia polyps and intraspecific similarity among Aurelia polyps. Therefore, our results suggest that the genotype of Aurelia polyps is an important factor in the structure of the symbiotic bacterial communities. Weiland-Bräuer et al. concluded that there were large differences in the microbial compositions of Aurelia polyps between North Sea/Roscoff and Baltic Sea subpopulations, similar to the present study (Weiland-Bräuer et al., 2015).

Animal hosts are known to be able to modulate their associated microbiome in terms of genotype regulation, such as via genetic expression or innate immune response activation (Borges, 2017). For example, in research using hydras as model animals, the variation in antimicrobial peptide genes in different species led to variation in the microbiomes among those species (Bosch, 2013). As part of the evolutionarily ancient marine phylum Cnidaria, Aurelia jellyfish may influence the structure of their microbiomes by producing various proteins or antimicrobial peptides via gene expression, which interferes with bacterial quorum sensing (selecting microbes from the environment or promoting the colonization and enrichment of targeted microbes) and inhibits bacterial colonization (Weiland-Bräuer et al., 2019). In addition, vertical transmission of bacterial communities from parents to offspring could be another reason for the differences between the bacterial communities of A. coerulea and A. solida polyps.

Hosts can conditionally shape the structure of microbial communities by interacting to maintain an overall stable state, especially cnidarians (Ainsworth et al., 2015). In the present study, the families Sphingomonadaceae (i.e., genera Sphingomonas), Rhizobiaceae (i.e., genera Phyllobacterium), and Burkholderiaceae (i.e., genera Ralstonia) were detected in all Aurelia polyp samples with high relative sequence abundances. Similarly, Weiland-Bräuer et al. proposed that the relative sequence abundance of Sphingomonadaceae was high in Aurelia collected in Roscoff, the North Sea, and the Baltic Sea throughout the life stages and decreased with the transition of Aurelia polyps to strobila to ephyra to juvenile medusae (Weiland-Bräuer et al., 2015). Peng et al. also detected a high abundance of Sphingomonadaceae and Burkholderiaceae in the associated bacterial community of wild Aurelia medusae (Peng et al., 2021). Several studies have suggested that cnidarian-associated bacterial communities are potentially involved in functional interactions and play a positive role in host-environment adaptation (Ziegler et al., 2017; Podell et al., 2020; Roach et al., 2020; Tong et al., 2020). Hence, the family Sphingomonadaceae may be a member of the core microbial communities associated with Aurelia populations worldwide. Further studies on the function of Sphingomonadaceae would help to gain insight into the impact of microorganisms on hosts.

Bacteria are functional components of most marine multicellular organisms, especially cnidarians (Morelan et al., 2019). Previous studies have agreed that the core microbiome was closely related to the health, growth, environmental adaptability, and production of the hosts (Bosch, 2013; Ainsworth et al., 2015; Brown et al., 2017; Weigel and Erwin, 2017; Weiland-Bräuer et al., 2020). In our study, a total of 20 core OTUs were identified in the bacterial communities of the scyphozoan body parts, predominantly the genera Sphingomonas, Ralstonia, and Phyllobacterium. The genus Sphingomonas (family Sphingomonadaceae), a dibenzofuran- and dibenzodioxin-degrading bacterium with potentially interesting properties for bioaugmentation of contaminated sites (Roggo et al., 2013), was the most abundant taxon in the microbiomes of all Aurelia polyps. Furthermore, Feng et al. (2017) identified that certain Sphingomonas spp. were potentially able to degrade chlorpyrifos, indicating that Sphingomonas could enhance the survival rate of Aurelia polyps in contaminated environments. The genus Phyllobacterium, which can produce exopolysaccharide (Li et al., 2017), was also dominant in Aurelia polyps. The exopolysaccharide could form protective barriers between cells and the environment, regulate cell growth and senescence, and affect cell division and differentiation (Flores-Felix et al., 2018), suggesting that the genus Phyllobacterium may be closely related to the transformation process between Aurelia life stages and host-environment adaptation.

Functional prediction is the first step in determining how microbiome biochemical processes affect the ecological functions of hosts (Ross et al., 2018). FAPROTAX is a promising tool for predicting ecologically relevant functions of bacterial and archaeal taxa derived from 16S rRNA amplicon sequencing (Louca et al., 2016). In this study, according to the FAPROTAX database, aerobic chemoheterotrophy in relation to C cycling was the primary function of Aurelia polyp-associated bacteria associated with numerous bacteria, such as Sphingomonadaceae and Microbacteriaceae. This demonstrated that Aurelia polyps were the main foundation of essential nutrients to support the microbial growth of the associated bacterial communities. Ross et al. (2018) reported that Aurelia medusae are potential bacterial vectors and may harm aquaculture activities, as their microbiomes harbor potential fish pathogens. Similarly, in this study, some potential animal pathogens (i.e., Coxiellaceae) or parasites (i.e., Bdellovibrionaceae and Haliangiaceae) were present in both A. coerulea and A. solida polyps. Moreover, nitrification and denitrification, two functions related to the N cycle that have been reported to be involved in waste removal in host coral (Beman et al., 2007), were potential functions that were abundant in both Aurelia polyp species. Hence, these bacteria may have similar N cycle-related functions in Aurelia polyps, assisting with the adaptation of polyps to nitrate stress in the ambient environment.

Each intraspecific group had unique bacterial taxa, which accounted for less than 1% of the relative sequence abundance. These unique bacteria are often overlooked because of their low abundance but may be critical to the functional maintenance of hosts (Shade et al., 2014). For example, Actinobacteria sp. and Ralstonia sp., located in zooxanthellae and coral intestinal epithelial cells, have low abundance and perform a vital role in the metabolism of the coral Acropora granulosa (Ainsworth et al., 2015). Unlike the relative stability of the core microbiome in Aurelia polyps ( Figure 1 ; Table 2 ), the unique microbes were less controlled and more sensitive to environmental variation. We speculated that unique microbes were preserved by the hosts from the native environment and potentially contributed to the environmental resilience of Aurelia polyps. Further study on unique microbes could comprehensively elucidate the mechanism of host selection for microbes in the environment.