Host Hybridization Dominates over Cohabitation in Affecting Gut Microbiota of Intrageneric Hybrid Takifugu Pufferfish

ABSTRACT Microbial symbionts are of great importance for macroscopic life, including fish, and both collectively comprise an integrated biological entity known as the holobiont. Yet little is known as to how the normal balance within the fish holobiont is maintained and how it responds to biotic and/or abiotic influences. Here, through amplicon profiling, the genealogical relationship between artificial F1 hybrid pufferfish with growth heterosis, produced from crossing female Takifugu obscurus with male Takifugu rubripes and its maternal halfsibling purebred, was well recapitulated by their gut microbial community similarities, indicating an evident parallelism between host phylogeny (hybridity) and microbiota relationships therein. Interestingly, modest yet significant fish growth promotion and gut microbiota alteration mediated by hybrid-purebred cohabitation were observed, in comparison with their respective monoculture cohorts that share common genetic makeups, implying a certain degree of environmental influences. Moreover, the underlying assemblage patterns of gut microbial communities were found associated with a trade-off between variable selection and dispersal limitation, which are plausibly driven by the augmented social interactions between hybrid and purebred cohabitants differing in behaviors. Results from this study not only can enrich, from a microbial perspective, the sophisticated understanding of complex and dynamic assemblage of the fish holobiont, but will also provide deeper insights into the ecophysiological factors imposed on the diversity-function relationships thereof. Our findings emphasize the intimate associations of gut microbiota in host genetics-environmental interactions and would have deeper practical implications for microbial contributions to optimize performance prediction and to improve the production of farmed fishes. IMPORTANCE Microbial symbionts are of great importance for macroscopic life, including fish, and yet little is known as to how the normal balance within the fish holobiont is maintained and how it responds to the biotic and/or abiotic influences. Through gut microbiota profiling, we show that host intrageneric hybridization and cohabitation can impose a strong disturbance upon pufferfish gut microbiota. Moreover, marked alterations in the composition and function of gut microbiota in both hybrid and purebred pufferfish cohabitants were observed, which are potentially correlated with different metabolic priorities and behaviors between host genealogy. These results can enrich, from a microbial perspective, the sophisticated understanding of the complex and dynamic assemblage of the fish holobiont and would have deeper practical implications for microbial contributions to optimize performance prediction and to improve farmed fish production.

ABSTRACT Microbial symbionts are of great importance for macroscopic life, including fish, and both collectively comprise an integrated biological entity known as the holobiont. Yet little is known as to how the normal balance within the fish holobiont is maintained and how it responds to biotic and/or abiotic influences. Here, through amplicon profiling, the genealogical relationship between artificial F1 hybrid pufferfish with growth heterosis, produced from crossing female Takifugu obscurus with male Takifugu rubripes and its maternal halfsibling purebred, was well recapitulated by their gut microbial community similarities, indicating an evident parallelism between host phylogeny (hybridity) and microbiota relationships therein. Interestingly, modest yet significant fish growth promotion and gut microbiota alteration mediated by hybrid-purebred cohabitation were observed, in comparison with their respective monoculture cohorts that share common genetic makeups, implying a certain degree of environmental influences. Moreover, the underlying assemblage patterns of gut microbial communities were found associated with a trade-off between variable selection and dispersal limitation, which are plausibly driven by the augmented social interactions between hybrid and purebred cohabitants differing in behaviors. Results from this study not only can enrich, from a microbial perspective, the sophisticated understanding of complex and dynamic assemblage of the fish holobiont, but will also provide deeper insights into the ecophysiological factors imposed on the diversity-function relationships thereof. Our findings emphasize the intimate associations of gut microbiota in host genetics-environmental interactions and would have deeper practical implications for microbial contributions to optimize performance prediction and to improve the production of farmed fishes. IMPORTANCE Microbial symbionts are of great importance for macroscopic life, including fish, and yet little is known as to how the normal balance within the fish holobiont is maintained and how it responds to the biotic and/or abiotic influences. Through gut microbiota profiling, we show that host intrageneric hybridization and cohabitation can impose a strong disturbance upon pufferfish gut microbiota. Moreover, marked alterations in the composition and function of gut microbiota in both hybrid and purebred pufferfish cohabitants were observed, which are potentially correlated with different metabolic priorities and behaviors between host genealogy. These results can enrich, from a microbial perspective, the sophisticated understanding of the complex and dynamic assemblage of the fish holobiont and would have deeper practical implications for microbial contributions to optimize performance prediction and to improve farmed fish production. function relationships of gut microbial communities thereof. Yet the mechanisms remain underexplored, and broad analyses of fish with variation in ecophysiology are particularly essential to obtain a generalized insight (26,27).
Takifugu (aka pufferfish) possesses a bioaccumulation capacity for lethal tetrodotoxin (28) and has a compacted genome replete with shortened intergenic and intronic sequences devoid of repetitive elements (29), rendering it a distinctive model organism. Captive Takifugu becomes nontoxic under domestication and is consumed as a delicacy widely across Asian countries (30)(31)(32). Among others, two sister species, namely, Takifugu rubripes (31) and Takifugu obscurus (33), have the most promising traits of interest, which are fast growth and hypotonic tolerance, respectively. Aquaculture production of both species has increased progressively to satisfy a growing market demand (34,35). Indeed, T. obscurus had been acclimated into complete freshwater residence (33), in contrast to T. rubripes that obligately resides in brackish water (36). Therefore, understanding how and determining the extent to which the growth performance can be improved and understanding how the environmental adaptability can be promoted have urgently become information needed to facilitate the rapid development of the pufferfish industry.
The gut microbiota contributes to the health and welfare of farmed fishes and can be regarded as an intermediate phenotype determined by both host genetics and environment, necessitating the incorporation of pairing the gut microbiota with host genetics for an assessment and prediction of trait performance (16,37). Previously, we showed structurally and functionally different gut microbiota among farmed juvenile T. obscurus at an inconsistent growth rate (38). Recently, we obtained viable artificial F1 hybrids produced by crossing T. obscurus ($) and T. rubripes (#) (39). Using a hybrid Takifugu pufferfish model, this study attempts to address the following questions. (i) Will the host genetics, environment, and interactions between them impact host growth and gut microbiota, and to what extent will these impacts occur? According to holobiont concepts (1-3), we hypothesized that changes in genetics and culture modes of the host will influence trait performance and/or microbial assemblage within the pufferfish holobiont. To address this hypothesis, the F1 hybrids together with their maternal halfsibling T. obscurus purebreds were reared separately by genealogy (monoculture) or equal proportionally cohabitated (coculture), and the host body mass and gut microbiota profile were analyzed integratedly. (ii) What are the phenotypic and functional implications of the gut microbiota underlying host genetics-environmental interactions? The propensity for colonization and proliferation in the gut varies greatly across microbial taxa (40,41), and as such, we hypothesized that the compositional and functional congruence of disturbed gut microbiota can be resolved. Accordingly, the bacterial phenotypes and predicted metagenomes of gut microbiota were inferred. (iii) What roles of different ecological processes govern the assemblage of the gut microbial community? Given the potential homogenizing effects (42) on the gut microbiota of cohabited fishes, we hypothesized that the interhost dissimilarity of the gut microbial community would be reduced while in cohabitation. To test this hypothesis, the patterns of ecological processes were discerned by the inference of phylogenetic signals and quantification of determinacy and stochasticity influence.

RESULTS
Synergistic effects of hybridization and cohabitation on fish growth: the host aspect of the pufferfish holobiont. The individual body weight and fork length were recorded to measure the explained variance in fish body mass by (i) host hybridity, i.e., hybrid ("H") versus purebred ("P"); and by (ii) culture modes, i.e., coculture ("co") versus monoculture ("mo") (see Fig. S1A and Table S1 in the supplemental material), respectively. Of note, such an unambiguous grouping of fish individuals by hybridity, particularly for the cocultures, is achieved by visually distinguishing the fish skin color patterns. That is, hybrids produced by T. obscurus ($) and T. rubripes (#) tend to be darker with more shattered yellowish-white spots spread over the dorsal part as described previously (39,43). First, we observed significantly larger body weight in hybrids than that of their maternal halfsibling purebreds in both "full" data sets with a sum of 400 pufferfishes (aligned rank transformed [ART] two-way analysis of variance [ANOVA], P , 2.2e-16) (Fig. 1A) and sub ones with 48 randomly selected individuals for downstream gut microbiota profiling (two-way ANOVA, P = 3.0e-12) (Fig. 1B). Similar trends held true for the fork length in both full data sets (ART two-way ANOVA, P , 2.2e-16) (Fig. 1C) and sub ones (two-way ANOVA, P = 6.4e-11) (Fig. 1D). These results indicated an evident growth vigor driven by host hybridization. In addition, significantly increased body weight and, to a lesser extent, fork length were observed in hybrid-purebred cohabitants compared with those of their respective monocultures that differed in genetic background for both full data sets ("weight," P = 5.5e-13; "length," P = 8.9e-6; ART two-way ANOVA) and sub ones (weight, P , 0.05; length, P = 0.066; two-way ANOVA), indicating an effect of cohabitation on pufferfish growth performance. Importantly, significant interactions between above two explanatory variables were also observed in both full data sets (length, P , 0.01; ART two-way ANOVA) and sub ones (weight, P , 0.01; two-way ANOVA), indicating a synergistic effect between host genetics and culture mode that affected fish growth.
Furthermore, an individual pufferfish was clustered on the basis of its aforementioned body mass measures by principal-component analysis (PCA) to obtain a combined metric of fish growth traits. Also, the quality of representation (cos 2 ) for both body mass measures was equally greater than 90% for the PC1 in either full data sets or sub ones. The per-individual PC1 (Table S1) was used as a body mass proxy for Body Length (cm) N=400 N=400

N=48 N=48
Genetics: P < 2e-16*** CultureMode: P < 5e-13*** Genetics ). Boxplots denote differences of fish body weight across groups within the full data set (A) and subset (B) and fish body length within the full data set (C) and subset (D). Statistics are shown as the panel titles. In each box, the horizontal bold bar denotes medians, the box height denotes interquartile range (25th to 75th percentile), and the whiskers denote the value range within 1.5 times the interquartile. Lowercase letters denote statistical significance reported at a confidence level of 0.95. Shown are two-dimensional scatterplots of all (E) and a subset (F) pufferfish individuals based on body mass metrics, including weight and length. The marginal density plots show the distributions of indicated body mass metrics per fish group. Lines denote the linear regression of pufferfish body mass (Spearman's rank).
Altered Gut Microbiota in Artificial Hybrid Pufferfish mSystems downstream analysis. To scrutinize and visualize individuals characterized by body mass, the results of univariate analyses of weight and length, respectively, were visualized by scatterplots, and all pufferfishes were well separated by both variables in full data sets (Fig. 1E) and sub ones (Fig. 1F).
Gut microbiota responses to host hybridization and cohabitation: the microbial aspect of pufferfish holobiont. Here, the microbiota of pufferfish gut and feed were profiled by 16S rRNA amplicon sequencing. A total of 6.4 million reads with an average length of 400 bp from 51 samples (a subset of 48 from 400 individuals; 3 technical feed replicates) were obtained. Through USEARCH-based read filtering and chimera discarding, followed by Ribosomal Database Project (RDP) classifier-based chloroplast sequence filtering, a total of 5,192 exact sequence variants (ESVs) were clustered. A sum of 5.8 million (90.7%) high-quality reads were retained (mappable to ESVs) with an average 113,561 reads per sample (Table S1). The coverage rate was estimated on rarefied data sets at 10,000 sequences per sample using rarefaction curves, which thereby ensured an adequate sequencing depth, and we retained all samples for further analysis ( Fig. S1B and C).
The calculated indices of metataxonomic diversity, including richness (observed ESVs), Chao1, Shannon, and Simpson (see Fig. S2 in the supplemental material), all suggested that the alpha diversity of the gut microbial community was markedly higher in purebreds (moP) than that in their maternal halfsibling hybrids (moH) in monoculture, and such a difference was slightly yet significantly diminished in hybrid-purebred cohabitants (coP versus coH) ( Fig. 2A). In addition, an anticorrelation between the alpha diversity and the fish body mass was observed (weight rho = 20.45, P , 0.01; length rho = 20.46, P , 0.001) (Fig. 2B). While, for each fish group separated by both genetics and culture modes, such an anticorrelation could be observed only in coculture purebred pufferfishes (length rho = 20.73, P , 0.01; weight rho = 20.69, P , 0.05) (Fig. 2B). Furthermore, comparisons of intercommunity dissimilarities (beta diversity) based on Jaccard distance showed that host body mass explained a larger fraction of variance (permutational multivariate analysis of variance [PERMANOVA], F = 4.66; R 2 = 0.089; P = 0.001), followed by host genetics (F = 1.65; R 2 = 0.031; P = 0.049) (Fig. 2C). Although the effect of culture mode alone was not statistically significant (F = 1.42; R 2 = 0.027; P = 0.123), its interactions with the host genetics was significant (F = 1.85; R 2 = 0.035; P = 0.03) (Fig. 2C). The permutational analysis of multivariate dispersion (PERMDISP) results showed that the intergroup dispersions between either genetics or culture modes were not significantly different in homogeneity (Permutest, P . 0.05) (see Fig. S3 in the supplemental material). Moreover, comparisons of weighted UniFrac distances showed that fish gut had greater dissimilarities to feed than those of all interhost pairs (P , 0.05; Mann-Whitney U test). Among others, the cocultured pufferfish of both genetics (coH and coP) had the shortest average intragroup distances and yet had the greatest dispersions (pairs "coH.coH" and "coP.coP"), suggesting a greater within-community similarity heterogeneity (Fig. 2D). Specifically, the gut microbiota was dramatically different compared with that of the feed, and it harbored less variation within genetics (asterisks in Fig. 2D) than between genetics (hashtags in Fig. 2D), which well recapitulated host genealogy. Furthermore, the discriminative taxa between fish groups were classified by the Random-Forest machine-learning approach, and a total of 23 bacterial families were identified to be able to classify a sample to group with 75% accuracy after 10-fold cross-validation (Fig. 2E). Wherein, two indicative bacterial families for hybrid pufferfish gut microbiota were identified, namely, Mycoplasmataceae (Tenericutes) and Fusobacteriaceae (Fusobacteria). Conversely, more different indicative families for purebred gut microbiota were characterized than the hybrid counterparts and brought into correspondence with the estimated higher alpha diversity of the formers. Consistently, the aforementioned two discriminative bacterial families were also found to be among the top-10 abundant taxa in hybrid pufferfish gut, by comparing the metataxonomic compositions of gut microbiota among fish groups ( Fig. 2F; see Fig. S4 in the supplemental material). Besides, a large majority of predominant gut bacterial families were found to be shared among fish groups, including Aeromonadaceae, Brevinemataceae, Flavobacteriaceae, and Vibrionaceae, in contrast to the predominant Enterobacteriaceae (48.9%), Moraxellaceae (24.9%) and Burkholderiaceae (22.2%) present in feed microbiota.
Phenotypic and functional implications of the gut microbiota underlying host genetics-environmental interactions. The organism-level phenotypes of bacterial communities, including Gram staining, oxygen tolerance, and pathogenic potential, were inferred by BugBase. The results showed that the gut microbiota of all fish groups are comprised of an overwhelming majority of Gram-negative bacteria (89%), which is significantly enriched in the purebred groups (95%) (Fig. 3A). Of note, the hybrid Altered Gut Microbiota in Artificial Hybrid Pufferfish mSystems aerobes (25%), and facultatively anaerobes (13%). While such proportions deviated slightly between fishes grouped by genetics, both anaerobes and facultative anaerobes were more present in the hybrid groups (42% and 19%, respectively) than those in the purebred ones (36% and 7%, respectively); albeit no statistical significance was found for the anaerobe phenotype (Fig. 3C). All fish groups were predicted to have prominent portions of bacteria that can form biofilms (71%), contain mobile elements (51%), and, to a lesser extent, be oxidative stress tolerant (44%) and potentially pathogenic (28%) (Fig. 3F to I).
To gain further biological insights into the functional differences among the gut bacterial communities that varied by fish groups, the metagenomes (KEGG pathway-level categories) were predicted by phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt). The results showed that both fish body weight and length were positively correlated with the abundance of the "carbohydrate digestion and absorption" pathway (Rho = 0.2 and 0.3) (Fig. S6). Moreover, certain pathways related to carbohydrate and lipid metabolisms were enriched in the gut microbiota of cohabitants (coH and coP) compared with their respective monocultured counterparts (Fig. S6). Specifically, two pathways, including "starch and sucrose metabolism" and "amino sugar and nucleotide sugar metabolism," with both being attributed to Brevinemataceae (Rho = 0.7 and 0.6), were enriched in the hybrid cohabitants (coH, log 10 linear discriminant analysis [LDA] of .3) (Fig. 3J). Whereas, other pathways, including "cysteine and methionine metabolism" and "other glycan degradation," which were attributed to Brevinemataceae (Rho = 0.6) and Rikenellaceae (Rho = 0.3), respectively, were enriched in the purebred cohabitants (coP, log 10 LDA of .2.5) (Fig. 3J). Regarding the monocultures, pathways, including "Nicotinate and nicotinamide metabolism" (metabolism of cofactors and vitamins) and "photosynthesis-antenna proteins" (energy metabolism), that were attributed to Brevinemataceae (Rho = 0.5) and some less abundant (average of ,1.5%) bacterial families (e.g., Rhodobacteraceae and Sphingomonadaceae; Rho, .0.8), respectively, were enriched in the monocultured hybrid ("moH," log 10 LDA of .2) (Fig. 3J). While, some pathways, including "transcription machinery" (genetic information processing) and "Arginine and proline metabolism," with both being attributed to Brevinemataceae (Rho = 0.5 and 0.6), were enriched in the monocultured purebreds (moP, log 10 LDA of .3) (Fig. 3J).
Ecological processes governing pufferfish gut microbiota assemblage in host genetics-environmental interactions. The phylogenetic signals of gut microbiota were determined by Blomberg's K statistic and Mantel's correlogram to quantify the underlying ecological processes governing microbial community assembly. A positive correlation between phylogenetic signal with a difference of fish body mass (PC1) was observed (P , 0.05), and yet it was merely within a narrow range of phylogenetic distance, rendering the quantification of phylogenetic turnover among closely related metataxa applicable to our data set. To detect ecological influences of various importance on submetacommunities that differed in relative abundance (0.01% cutoff), the gut microbiota ("all," 4,908 ESVs) were fractionated into "abundant" (.0.01%, 326 ESVs) and "rare' (,0.01%, 4,582 ESVs) submetacommunities. The calculated unweighted-standardized effect size of the mean nearest taxon distances (ses.MNTD) were well below zero for each of these submetacommunities per fish group (P , 0.01) (see Fig. S7 in the supplemental material), indicating the importance of environmental filtering in driving phylogenetic clustering of the gut metacommunity. Notably, the abundant submetacommunity in hybrid groups had a greater absolute unweighted ses.MNTD than that in the purebred counterparts (Fig. 4A), and such differences were anticorrelated with the host body mass (Fig. 4B), suggesting a lower degree of phylogenetic relatedness in the gut microbiota of the slow-growth purebred.
Furthermore, the pairwise phylogenetic beta-diversity (bMNTD) was estimated and compared against the null model to calculate the unweighted-beta nearest taxon index, bNTI, of which the absolute values greater than 2 denote significant phylogenetic turnover comparing with those of the randomized community. In addition, results of the Mantel's test showed significant correlations between bNTI and changes of body mass for each of the submetacommunities grouped by abundance, particularly for the abundant ones ( Fig. 4C; see Fig. S8 in the supplemental material). These results demonstrated that the gut microbiota assemblage was attributed to the shift from a small fraction of homogeneous selection (1;14%) to a great majority of variable selection (45 to ;61%) as the changes of host body mass increased, regardless of abundance (Fig. 4D, top). Regarding the stochastic ecological processes with an absolute bNTI value less than 2, their pairwise distances were further partitioned by capping absolute values of both bNTI and RC bc at 2 and 0.95, respectively. The results showed that dispersal limitation (19 to 31%) and "undominated" (3 to 26%) processes contributed largely to nondeterministic ecological processes for each subcommunity grouped   a a a a a   a a a a a a a   a a a a a a a   a a a a a  The top three stacked bars depict the relative importance of different ecological processes in governing microbial community assemblage for "all," "rare," and "abundant" subcommunities, respectively. And the bottom four stacked bars depict the "abundant" subcommunities of each indicated pufferfish group. Percent font colors, white and black, within the stacked bars denote determinacy and stochasticity fractions, respectively. All statistics were shown as subpanel titles.
Altered Gut Microbiota in Artificial Hybrid Pufferfish mSystems by abundance, while the homogenizing dispersal had little to no contribution (0 to 2%) (upper panel in Fig. 4D). Given that the phylogenetic signals of abundant subcommunities were significantly different between purebreds and hybrids (Fig. 4A), the influence of ecological processes determining the assembly of abundant subcommunities was further estimated per fish group. The results showed that the ratio between deterministic and stochastic processes was near six times larger in monoculture purebreds (76%/24%) than that in monoculture hybrids (35%/65%) (Fig. 4D, bottom). Moreover, such differences diminished slightly in cocultures, as evidenced by the decreased and increased proportions of variable selection in monoculture purebreds (from 76% to 62%) and monoculture hybrids (from 29% to 48%), respectively. Consistently, neutral modeling indicated greater deviations in purebred groups ('coP': R 2 = 0.61; 'moP': R 2 = 0.48) than those in hybrid groups ('coH': R 2 = 0.70; 'moH': R 2 = 0.69) while considering all metataxa (see Fig. S9 in the supplemental material).

DISCUSSION
The main objective of this study was to explore the impacts of host genetics, environment, and interactions between them on pufferfish growth and gut microbiota. Here, the gut microbiota of artificial F1 hybrid pufferfish with marked growth heterosis, along with their maternal halfsibling purebreds, were profiled by 16S rRNA amplicon sequencing, and the underlying interrelations of gut microbiota, host body mass, and fish culture modes were determined. Our results showed that the gut microbiota are compositionally and functionally distinct by host genealogy. Such differences of host growth and gut microbiota both were slightly and yet significantly altered in hybridpurebred cohabitants, compared with respective monoculture cohorts that shared common genetic make-ups. Furthermore, the phenotypic and functional differences and the ecological processes governing community assemblage of gut microbiota were inferred. Our results not only can enrich, from a microbial perspective, the sophisticated understanding of the complex and dynamic assemblage of fish holobionts but also will provide deeper practical implications for microbial contributions for optimizing performance prediction and improving the production of farmed fishes.
Effects of hybridization on pufferfish holobiont: implications of host metabolic priority shift toward carbohydrate utilization. Often, interspecific hybridization can give rise to progeny with enhanced or diminished viability and/or fertility, thereby conferring fitness superiority or inferiority, respectively, relative to parents (19,44). Here, the F1 hybrid pufferfishes exhibited superior growth performance and intact hypotonic acclimation capacity relative to their maternal halfsibling purebreds (Fig. 1), indicating enhanced fitness associated with intrageneric hybridity. Furthermore, this hybrid cohort can produce viable F2 progeny by backcrossing with the maternal parents equally well as the purebreds (data not shown), suggesting an unimpaired fertility. Considering the husbandry condition, pond management and feeding regime all were identical (within cocultures) or nearly so (between monocultures) for the fish tested here, and the above results indicated strongly that the intrageneric hybridity, over the environment if not entirely, confers growth vigor without apparent fertility cost. In fact, natural hybridization between a given pair of Takifugu spp. in the estuarine and near-shore habitats had been reported previously (30,45,46), demonstrating an ecologically and behaviorally incomplete prezygotic reproductive isolation thereof. Importantly, Takifugu spp. exhibit a relatively lower degree of intrageneric genetic differentiation (i.e., less genetically incompatible), as is evident by its well-conserved karyotype and synteny (46)(47)(48), rendering possible the postzygotic synapsis, recombination, and segregation of chromosomes in hybrid progenies. Moreover, artificial hybridization in Takifugu had been proven feasible and amendable (39,43), justifying its promising application potential in breeding practice. Although population genetics had revealed that clear genomic divergence correlated with hypotonic adaptation of T. obscurus (33) and fast growth of T. rubripes (31), compared with other Takifugu spp. with respective lower trait values, the underlying genetics mechanisms of heterosis warrant finer-scale investigations. One plausible explanation for such growth vigor is Altered Gut Microbiota in Artificial Hybrid Pufferfish mSystems paternal effects, as is evident by the comparably higher levels of growth hormones in both hybrid progeny (T. obscurus $ x T. rubripes #) and male parent (T. rubripes) than those in the female parent (T. obscurus) (43). Nevertheless, it remains elusive whether such heterosis is conferred by the allele-specific overdominance and/or by the genome-wide heterozygosity per se.
In addition to the extrinsic environmental factors, the hybrid fitness can be determined by the compatibilities and interactions, not only between the allospecific host DNA that have never before coadapted and coevolved together (discussed above) but also between their respective symbiotic microbiota that assembled and succeeded independently (49)(50)(51). Previously, we showed the T. obscurus gut microbiota was replete with copiotrophic bacteria (with higher numbers of rRNA genomic copies) that can proliferate faster, in contrast to the oligotrophic ones predominated in pond water and sediment, indicating a substantial niche barrier between pufferfish gut (nutrient rich) and the surrounding environment (nutrient poor) (38). In addition, three bacterial families, namely, Brevinemataceae, Rikenellaceae, and Mycoplasmataceae, were characterized to be autochthonous (core) pufferfish gut microbiota (38). Consistently, these three bacterial taxa, albeit varied in abundance, were also present in the gut microbiota tested here, regardless of host genealogy (Fig. 2F). Moreover, the hybrids did have a markedly different gut microbiota, as evidenced by a lower alpha diversity (species richness, Fig. 2A and B), higher beta diversity (community dissimilarity, Fig. 2C and D), disproportionate bacterial phenotypes and functional compositions (Fig. 3A, B, E, and J), higher phylogenetic relatedness (ses.MNTD, Fig. 4A and B), and more stochasticity influenced (Fig. 4D) than the purebreds. These results together with our previous findings indicated an evident host genetics-associated niche effect on the gut microbiota, which can be scaled up from the host genus (38) down to the species level (this study). Such a linkage between host genetics and gut microbiota is, albeit to a limited extent, reminiscent of the phylosymbiosis concept, which suggests that more phylogenetically related hosts tend to have more compositionally similar microbiota (8). Mounting evidence has emphasized the importance of both vertical and horizontal transmissions of microbiota in shaping symbiosis, which are underpinned by microbial acquisition and cultivation on or in the host (7,8). Considering a shared maternal source and hatching environment of the fish tested here, and almost negligible albeit interspecifically different contribution of paternal gamete in early life microbiota exposure thereof, the common microbiota between the hybrid and purebred may signify vertical transmission from the shared maternal parents. While such parent-to-progeny transmission is highly unlikely as faithful as seen in mammals (7,8), given the proneness to environmental disturbance of fish external fertilization. Conversely, fish gut microbiota can be recurrently sourced from the environment by horizontal transmission, which is determined by the gut physiochemical condition defined by host genetics (21) and/or by the microbial propensity for colonization and proliferation in a given gut environment (40,41).
Here, the pufferfish gut niches were generally characterized to be favored by the Gram-negative, anaerobic, and biofilm-forming bacteria, wherein the hybrid groups appeared to be more favored by the Gram-positive and facultatively anaerobic bacteria than the purebred ones, as indicated by Bugbase inference (Fig. 3A to I). Moreover, several carbohydrate metabolism-related pathways of the predicted gut metagenomes were enriched in the hybrid groups, which contrasted with the enriched lipid and amino acid metabolism-related pathways in the purebred groups (Fig. 3J). Previous studies revealed higher growth hormone levels in hybrid pufferfish (T. obscurus $ x T. rubripes #) than those in purebreds (T. obscurus), which can accelerate feed intake and feed efficiency and impact digestive metabolism of feed intake (43). In fact, growth hormone transgenesis in diverse fish species can cause a shift in host metabolic priority, wherein the dietary carbohydrate was used as a preferential energy source over protein and lipid, and its carbohydrate utilization was elevated to meet the high energy demand and spared for fast growth-associated biosynthesis of the protein and lipid (52)(53)(54). As such, we speculated that the observed metataxonomic and functional differences of gut microbiota likely reflect the disparate metabolic priorities between pufferfish hosts that differed in genealogy. Considering a likely increased feed intake and energy demand of the hybrid pufferfish tested here, more dietary carbohydrate and hydrolysis intermediates in the gut digesta could be expected, thereby fueling the growth of restricted bacterial taxa specialized for carbohydrate fermentation, such as Fusobacteriaceae and Mycoplasmataceae (i.e., taxa with high abundance). This point can be further strengthened by the positive correlation between the abundance of the "carbohydrate digestion and absorption" pathway and fish body mass ( Fig. 3J; Fig. S6), which can be observed more obviously in hybrid pufferfishes. Conversely, in the absence of a high growth hormone level as in the hybrid, the basic catabolic processes of dietary protein and lipid may have been maintained in the purebred pufferfish gut, from which some catabolites, including essential amino and lipid acids and derivatives, can be made available to and be extensively used by a wider range of microbes as growth factors or precursors, thereby nourishing growth and increasing the overall alpha diversity thereof. In turn, the different gut microbiota of hybrid and purebred groups may also alter the gut physiochemical environment, for instance by microbial metabolites, and further enlarge such intercommunity dissimilarity. Altogether, the above results demonstrated a profound influence of artificial intrageneric hybridization on the pufferfish holobiont, which is plausibly associated with a shift in the metabolic priority of the host.
Influences of cohabitation on the pufferfish holobiont: implications of environmental enrichment for fish. Practically, by growing different fish species or populations together, polyculture or coculture not only can enrich ecosystem biodiversity and provisioned function (55) but also enhance environmental enrichment by providing mechanical and physiological stimuli to the cultured fish (56). Here, coculturing hybrid and purebred pufferfishes caused changes to both fish growth (Fig. 1) and gut microbiota (Fig. 2), compared with the equally dense monoculture counterparts with the same genetic background. Despite the dispersions due to an insufficient sample size, the increasing tendency of fish growth in the cocultured subset (sub) remained robust, particularly in the purebred group (Fig. 1). By increasing biodiversity, for instance, by complementarily integrating multitrophic or monotrophic aquaculture, polyculture, or coculture can promote fish growth and reduce waste discharge, ultimately expediting the transition toward more environmentally sustainable aquaculture (57). Nevertheless, the fish growth promotion by coculture observed here is less likely as a result of trophic complementation (e.g., feed utilization) since both parental species, i.e., T. rubripes and T. obscurus, are omnivores (58). Few promising behavioral studies had shown that T. rubripes exhibits quicker acclimation in new environment and shorter latency to both feed and novel object, suggesting a distinct proactive behavior that may linked to fast-growth (59). The genetic basis underpinning such proactiveness was found to be controlled by a combination of coadapted multiallelic loci together with small pleiotropic effect (59). Such an intrageneric difference of host behavior would be a likely explanation for the growth enhancement by cohabitation observed here. It is tempting to speculate that, as a way of environmental enrichment, coculturing of hybrid and purebred pufferfishes that differed in stress copying styles or personalities can sufficiently intensify social interactions (56). For instance, by providing sensory stimuli and social novelty and meanwhile maintaining conspecific and/or congeneric familiarity, social interactions can stabilize the social hierarchy of the pufferfish cohabitant and mitigate the stress response that will inhibit activity and feeding, thereby enhancing growth and ultimately optimizing fitness and survival, as is seen in other fish coculture systems (57).
According to the metacommunity theory, the host-associated microbiota does not exist in isolation and permanence but rather can change membership by reciprocally transmitting in and out from other hosts and the environment (7,8). Here, the purebred-hybrid cohabitants both have significantly less dissimilar microbiota either between or within genetic groups than that of respective monocultures (Fig. 2D), indicating a measurable influence of cohabitation by diminishing interhost dissimilarity. Of note, an inflated dissimilarity dispersion with reduced distance in each community pair, including cohabitants, was also observed (longer interquartile range in Fig. 2D). Such results somehow reflect the augmented interhost heterogeneity of gut microbiota by cohabitation, in a tendency of reducing dissimilarity particularly apparent among cohabitants (i.e., rows coH.coP, coH.coH, and coP.coP). Such homogenizing effects appeared to have a subtle effect on the phylogenetic relatedness of gut microbiota, as indicated by the comparable ses.MNTD indexes between culture modes per genetic group (Fig. 4A). These results suggest that cohabitation does not homogenize microbial communities simply by averaging diversities (richness and similarity) but rather it involves complex interactions within a metacommunity across all hosts and its environment. As aforementioned, the two Takifugu spp. discussed here vary considerably in behaviors (59), and this variation may have implications not only for host growth promotion (e.g., by enriched social interactions) but also for community membership changes of microbes (e.g., by in-between transmission). By stimulating host behavioral and physiological changes, animal social interactions had been shown to have multitude effects on the assemblage of symbiotic microbiota (7,56,60). According to the social microbiome concept (7), animal microbiota not only rely on the host for resources but also depend on social interactions and networks of the host (particularly group-living animals) to be transmitted. Therefore, proximity and connectivity within a common environment shared by cohabitated animals (in our case hybrid-purebred cohabitants) are expected to probabilistically promote direct (physical) or indirect microbial transmissions (e.g., fecal contamination) and dispersal, ultimately reducing community dissimilarity (7). Importantly, indirect transmission may involve transient environmental persistence that selects for microbes capable of survival as free-living organisms, for instance, conferred by oxygen tolerance, motility, and surface colonization (7,41). In supporting this notion, the two discriminative taxa Fusobacteriaceae and Mycoplasmataceae (Fig. 4E and F) in hybrids by either culture modes both have a strict or facultative anaerobic lifestyle and fermentative abilities, indicating a lack of extrahost persistence as seen in other fish species (61). Conversely, among the most prevalent taxa shared by all pufferfishes, Aeromonadaceae, Flavobacteriaceae, and Vibrionaceae (Fig. 3F) are characterized to be aerobic, motile, and surface adhesive, which is consistent with their ubiquity in the aqueous environment and active interhost transmissions (62,63).
Furthermore, we observed distinct assemblage patterns between the gut microbiota of purebred and hybrid groups within coculture. For instance, the influence of determinacy in the monocultured purebred group was weakened in coculture (from 76% to 62%), which contrasted with that of the hybrid counterparts that were enhanced in coculture compared with monoculture (from 35% to 48%) (Fig. 4D), suggesting slightly converged ecological processes governing community turnover for both hybrid and purebred cohabitants. Among the determinacy, variable selection (bNTI, .2), rather than homogeneous selection (bNTI, ,22), varied from moderate to strong as the body mass differences increased (Fig. 4C), highlighting the importance of augmented spatial environmental heterogeneity in determining gut microbial community assemblage (6,10). In our case of coculture, such variable selection could be attributed to the enhanced complexities of host-microbial interactions within and/or between hybrid and purebred groups. Given the cohabitants with presumptively different behavior and the spatially and temporally connected gut microbiota, either or both can cause more fluctuations in the gut microenvironment, thereby imposing an enlarged heterogeneous influence on the community composition, as reflected in increased dissimilarity dispersion with reduced dissimilarity within coculture community pairs (Fig. 2D). Moreover, a large fraction of stochastic influence was found governed by dispersal limitation (;30%) and undominated processes (;15%) and little (,1%) to no influence by homogenizing dispersal (Fig. 4D). The undominated fraction can be regarded as consequences of attenuated selection and/or dispersal, with both being the opposite of variable selection (6). In supporting this notion, the undominated fraction is absent only in the gut microbiota of the monocultured purebred group, which is influenced largely by variable selection. Notably, the stochastic influence, particularly dispersal limitation, was found prominent in monoculture hybrids (44%) but less so in cocultured counterparts (27%). In contrast with homogenizing dispersal (stronger dispersal), dispersal limitation indicates the host effect (e.g., as patched niche), by which the magnitude to which microbes transmit among individual gut was constrained (6). Such higher levels of stochastic influence in hybrids, particularly monocultures, may be associated with a lower degree of variations in intra-and/or extrahost environments, which were implicated with prioritized metabolism of carbohydrate and a lack of environmental enrichment for the host (discussed above). In line with this notion, a neutral model-based inference suggested a greater importance of passive migration and demographic processes in explaining the compositional turnover for both cocultured and monocultured hybrid groups, as reflected by the higher model fitness (R 2 , ;0.7) and yet slightly lower migration rate (m, ;0.05), compared with that of purebred ones (R 2 , ;0.5; m, ;0.18). These results stress that weaker dispersal (i.e., higher fractions of dispersal limitation) can cause dissipation within local gut communities in the form of reducing alpha diversity (9), as seen in the hybrid group ( Fig. 2A).
In conclusion, we show that the gut microbiota relationships between artificial F1 hybrid and purebred pufferfishes well recapitulated their genealogy, indicating evident host genetic effects on gut microbiota. Furthermore, significant fish growth promotion and gut microbiota alteration mediated by hybrid-purebred cohabitation were observed, in addition to the influence of intrageneric hybridization to a larger extent. The compositional and functional congruence of gut microbiota between host genealogy were likely determined in part by the host metabolic priority shift toward carbohydrate utilization. Moreover, the underlying assemblage mechanisms of gut microbiota were found involved in a trade-off between variable selection and dispersal limitation, which was plausibly driven by augmented social interactions between hybrid and purebred cohabitants with different behaviors. Nevertheless, this study cannot explicitly ascertain whether such hostgenetic determination is inheritable in the long run and to what extent and if the hybrid heterosis can remain intergenerationally constant (i.e., without breakdown), considering the lack of parental samples, particularly of brackish water-resident T. rubripes, examined in parallel. An important caveat in this study is the lack of replicated culture ponds due to limited space and resources. Future studies with enlarged pond replicates and sample sizes need to account for genetic and ecophysiological aspects of the holobiont that contribute to the host fitness and microbial transmission. Altogether, our findings emphasized the critical interfacial roles of gut microbiota in host genetics-environmental interactions and would have deeper practical implications for microbial contributions to optimize performance prediction and improve the production of farmed fishes.

MATERIALS AND METHODS
Experimental design and sampling. The 180 days posthatch (dph) juvenile pufferfishes were obtained from a single commercial supplier (Jiangsu Zhongyang Group Co., Ltd.) located in Nantong, East China, in September 2020, as described previously (39). In brief, the artificial F1 hybrids (H) were produced by crossing four each of T. obscurus ($ siblings) and T. rubripes (# siblings), and their halfsibling T. obscurus purebreds (P) were produced by the same maternal parents. The hybrids and purebreds were either (i) grouped by the same genealogy and then reared separately (2Â monoculture, moH and moP) or (ii) equal proportionally cohabitated (1Â coculture, coH and coP), in a total of three adjacent indoor ponds (6 m by 6 m by 1 m) at an equal stocking density (Fig. S1A). The husbandry conditions, pond management, and feeding regime all were controlled to be identical (within cocultures) or nearly so (between monocultures). A sum of 400 individual pufferfishes, more precisely 100 fishes per group, were captured randomly from each pond and were immediately weighted and measured for fork length. For gut microbiota profiling, a subset of 48 individuals (n = 12 per groups) from aforementioned fishes were selected randomly and sacrificed. After being surface sterilized using 75% ethanol, the distal intestine (2 cm long) was removed with a sterilized spatula to collect the gut content. Subsequently, all gut content samples were placed in separate sterile tubes, and the feed samples were packaged in airtight sterile plastic bags. All microbial samples were preserved with dry ice during shipment and stored at 280°C. The detailed metadata for each sample were recorded in Table S1. The sampling was performed in line with the protocol approved by the Ethics Committee of Experimental Animals (Hohai University) and was in accordance with the Animal Care Guidelines issued by the Ministry of Science and Technology (China).
Altered Gut Microbiota in Artificial Hybrid Pufferfish mSystems DNA isolation and 16S rRNA amplicon sequencing. The microbiota of gut contents (n = 48) and feed (n = 3) were profiled by 16S rRNA amplicon sequencing. In brief, DNA was extracted with the FastDNA spin kit for soil (MP Biomedicals, Irvine, CA) and verified with 0.8% agarose gel electrophoresis. DNA integrity and concentration were evaluated by a Nanodrop ND2000 instrument (Thermo Fisher Scientific Inc., Waltham, MA) and a Quant-iT PicoGreen double-stranded DNA (dsDNA) kit (Thermo Fisher Scientific Inc.), respectively. The V4/V5 hypervariable region was amplified using primer pair 515F (59-CTGCCAGCMGCCGCGGTAA-39) and 926R (59-CCGTCAATTCMTTTRAGTTT-39) barcoded with unique eightnucleotide sequences per sample. The 25-mL PCR system included 40 ng template DNA, 5 mL 5Â reaction buffer, 5 mL 5Â GC buffer, 2 mL dNTP (2.5 mM), 10 pM barcoded primer pairs, and 0.75 U Q 5 high-fidelity DNA polymerase (New England BioLabs [NEB], Ipswich, MA). All samples were tested in triplicate together with no-template controls. The thermocycling conditions were set at 98°C for 2 min; 25 cycles at 95°C for 15 s, 55°C for 30 s, and 72°C for 30 s; followed by a final step of 72°C for 5 min. Amplicons were subjected to electrophoresis with a 2% agarose gel, were purified by the AMPure XP kit (Beckman Coulter GmbH, USA), and were quality checked by a high-sensitivity DNA kit (Agilent, Santa Carla, CA) and Quant-iT PicoGreen dsDNA kit. Purified amplicon DNA libraries were constructed using a TruSeq Nano DNA LT library prep kit (Illumina, San Diego, CA), and the 400-bp insertion was paired-end sequenced with the MiSeq reagent kit V3 chemistry on an Illumina MiSeq platform.
Read processing and microbiota profiling. Raw amplicon reads were trimmed by Trim-Galore (Cutadapt v.1.4.2 and FastQC v.0.10.1) with a default Q20 quality cutoff (64). Trimmed reads were processed according to a pipeline described previously (65) with USEARCH v.10.0 (66). Briefly, paired-end reads were merged by "-fastq_mergepairs"; primer sequences were removed by "-fastx_truncate"; and low-quality and redundant reads (error rates, .0.01) were filtered and dereplicated by "-fastq_filter" and "-fastx_uniques," respectively. The ESVs were clustered after filtering chimeras by "-unoise3." Subsequently, the ESV table was generated by "-otutab," and taxonomically assigned by the sintax algorithm with the Ribosomal Database Project (RDP) classifier (RDP training set v.16) (67). The potential organelle (mitochondria and chloroplast)derived reads were removed by taxonomy-based filtering. The sequencing depth was normalized by "otu-tab_norm," yielding 10,000 sequences per sample. Community diversity was profiled by "-alpha_div," "-alpha_div_rare," "-cluster_agg," and "-beta_div." To identify the featured taxa among different microbial communities, RandomForest (v.4.6-14) (68) was used to classify the microbial taxa in the family level across different groups with parameters of "ntree = 1000, importance = TRUE, proximity = TRUE," followed by cross-validation by the rfcv() function. For functional inference, the ESV table was normalized according to 16S rRNA copy numbers and annotated by PICRUSt (69) with closed-reference picking against a compatible version of the Greengenes database (v.13_5) (70). Differences in relative KEGG pathway abundance were identified with a cutoff of logarithmic LDA score greater than 2 combined with the linear discriminant analysis effect size (LEfSe) algorithm (71). Bacterial phenotypes, such as oxygen tolerance, Gram staining, and pathogenic potential, within each indicated microbial community were inferred by BugBase (72) with the taxonomic level set at bacterial family.
Quantifying ecological processes for governing microbial community assemblage. The tree distance-based phylogenetic relatedness inference of a given group of gut microbial community was calculated by phyloSignal() and phyloCorrelogram() within the "Phylosignal" package (73). The standardized effect size of the mean nearest taxon distances (ses.MNTD) was calculated by the ses.mntd() function within the "Picante" package (74) per community, with options of "taxa.labels" as the null expected model, unweighted, and 999 randomization. The ses.MNTD that were significant less than zero indicate closer phylogenetic relatedness than expected by chance (75). To further quantify the relative importance of ecological processes in governing assemblage in microbial community, both metrics of b-nearest taxon index (bNTI) and Bray-Curtis-based Raup-Crick (RC BC ) were calculated to evaluate the contributions of determinacy and stochasticity, respectively (6). In doing so, the phylogenetic beta diversity was calculated based on unweighted intercommunity MNTD metric (bMNTD) by the comdistnt() function. Thereafter, the deviation of observed bMNTD from the null expectation (999 randomization), i.e., bNTI, was computed as described previously (6). The influence of determinacy, including homogeneous selection and variable selection, was inferred with significance bNTI that less than 22 and greater than 2 (6). Subsequently, the nonsignificant fractions that have an jbNTIj value less than 2, reflective of stochasticity, were further subdivided according to b-diversity metrics of RC BC , wherein the RC BC indexes normalize the deviation of observed Bray-Curtis from null expectation (9,999 randomization) into a range between 21 and 1. The significant RC BC indexes that are less than 20.95 or greater than 0.95 were interpreted as influenced by homogenizing dispersal and dispersal limitation, respectively. The nonsignificant fractions (i.e., jRC BC j of ,0.95) were considered community assemblage influenced by undominated processes without a main effect by any given ecological pressure (i.e., weak dispersal and/or selection).
Statistical analysis. All statistical analyses were conducted under the R environment (v.3.6.0). Prior to the analysis, the normality of residuals and homogeneity of variances were tested by shapiro.test() and bartlett.test(), respectively, followed by parametric or nonparametric tests. In doing so, the group means of parametric data sets were compared by one-way and/or two-way (i.e., plus asterisk between two explanatory variables in the formula to account for potential interactions) analysis of variance (ANOVA) using aov(), followed by post hoc test by TukeyHSD(). Alternatively, the rank sums of nonparametric data sets were compared by Kruskal-Wallis test with kruskal.test(), followed by Wilcoxon rank-sum test by wilcox_test(). The two-way nonparametric statistical tests were performed after the ART procedure by art() within the "ARTool" package, followed by model fit test by anova(). Statistical significance was considered while reported false discovery rate (FDR)-adjusted P values were less than 0.05. The PERMDISP test was used to evaluate the homogeneity of multivariate dispersion of groups by betadisper (), followed by permutation test with permutest(). Principle-coordinate analyses (PCoAs) were conducted based on the Jaccard dissimilarity metric by cmdscale(), followed by the permutational multivariate analyses of variance (PERMANOVA) test on the model considering fish body mass (PC1) as a covariate plus the interactions between host genetics and culture mode by adonis() within "vegan" package (76). The clustering analysis on fish body mass was visualized after principal-component analysis by prcomp(). Potential correlations between bNTI and indicated trait difference (Euclidean-based trait distance matrix) were assessed by mantel() within the "ecodist" package (77) with 9,999 permutations. The fitness of the neutral model (78) on the gut microbial community was evaluated by neutral.fit() within the "MicEco" package.
Data availability. All raw reads generated from this study were deposited at the National Center for Biotechnology Information (NCBI) sequencing reads archive (SRA), under BioProject accession number PRJNA841699.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.