Microbial embryonal colonization during pipefish male pregnancy

While originally acquired from the environment, a fraction of the microbiota is transferred from parents to offspring. The immune system shapes the microbial colonization, while commensal microbes may boost host immune defences. Parental transfer of microbes in viviparous animals remains ambiguous, as the two transfer routes (transovarial vs. pregnancy) are intermingled within the maternal body. Pipefishes and seahorses (syngnathids) are ideally suited to disentangle transovarial microbial transfer from a contribution during pregnancy due to their maternal egg production and their unique male pregnancy. We assessed the persistency and the changes in the microbial communities of the maternal and paternal reproductive tracts over proceeding male pregnancy by sequencing microbial 16S rRNA genes of swabs from maternal gonads and brood pouches of non-pregnant and pregnant fathers. Applying parental immunological activation with heat-killed bacteria, we evaluated the impact of parental immunological status on microbial development. Our data indicate that maternal gonads and paternal brood pouches harbor distinct microbial communities, which could affect embryonal development in a sex-specific manner. Upon activation of the immune system, a shift of the microbial community was observed. The activation of the immune system induced the expansion of microbiota richness during late pregnancy, which corresponds to the time point of larval mouth opening, when initial microbial colonization must take place.

Microbes can be parasites that have a detrimental impact on their hosts. However, a vast diversity of non-pathogenic microbes lives as commensals or mutualists on their hosts. Altogether they constitute the "holobiont" [1][2][3][4] . As specialized entities within organs 5 , microbiota shape almost every aspect of animal physiology (behaviour, reproduction, fitness) and may even play a role in hybridization and speciation [6][7][8][9][10] . For the successful host development, microbiota are indispensable, as they foster the production of polysaccharides and vitamins, boost the maturation of the host immune system, protect the host against invasions by pathogens, and maintain host tissue homeostasis [11][12][13][14][15][16] .
In contrast to maternal and paternal genes that are inherited in a Mendelian manner, the microbiome was originally acquired from the environment. Its transfer from parents to offspring revitalizes aspects of Lamarckism 4 . Siblings have a more similar microbiome than unrelated newborns, and a litter harbours offspring with more similar microbiota than offspring from different litters [17][18][19][20] . This implies that both the environment and the genetic background shape the early microbiome composition 21,22 . The maternal environmental experience may thus modulate the embryonic immune system in two ways: via the direct transfer of immunity 23 and via the transfer of commensal bacteria 24 . This suggests a co-adaptation between parental investment, the host microbiota and the immune system.
Several features make the sex-role reversed pipefish Syngnathus typhle ideal for studying microbial transfer from parents to offspring. While pregnancy has evolved multiple times independently in most vertebrate groups 25,26 , syngnathids (pipefishes and seahorses) are the unique representative of male pregnancy evolution 27,28 . Mothers transfer their eggs into the paternal brood pouch, which is connected to a placenta-like system, where they are bred for 4-6 weeks 28,29 . We predict that in pipefish both vertical maternal transmission of specific commensal bacteria and additional paternal transmission during male pregnancy may shape initial translocation of embryonic microbiota and influence immunological maturation. This system permits the disentangling between transovarial microbial transfer (maternal) and microbial contribution during pregnancy (paternal), two traits that are in all other viviparous systems intermingled within the female.
SCiEntiFiC RePORTS | (2019) 9:3 | DOI:10.1038/s41598-018-37026-3 microbiota differed from the paternal brood pouch microbiota ( Fig. 2A). Brood pouch microbiota of non-pregnant males differed from both the gonad and pregnancy associated microbiota, while it was more similar to non-pregnant brood pouch microbiota ( Fig. 2A). Microbiota in the brood pouch of late-pregnant males was distinct from earlier pregnancy stages ( Fig. 2A). In addition, brood pouch microbiota during late stage pregnancy showed a distinct clustering of microbiota, when their parents had experienced an injection with heat-killed bacteria (Fig. 2B). The most abundant Otu00001: Oceanospirillaceae, Marinomonas with 28% was overrepresented in late pregnancy brood pouch microbiota ( Fig. 2A), in particular if fathers were exposed to prior immune challenges (Fig. 2B). In contrast, the second and third most abundant OTUs Halomonadaceae, Halomonas with 10.3% and Bacillaceae_1, Aeribacillus with 6.2% were predominantly prevalent in the gonads of injected females (GO_I) (Fig. 2B). The microbiota in the brood pouch of non-pregnant, early, mid and late pregnant males were, independent of the paternal immunological activation, more uniform than the microbiota in the female gonads (Fig. 2). This could imply that eggs transferred from mothers into the paternal brood pouch came in contact with a novel microbial community, which is during late pregnancy subject to drastic changes.
Bacterial α -diversity within and between developmental stages and the impact of parental immunological activation. For the assessment of bacterial community shifts within and between developmental stages on the surface of pipefish eggs/embryos during male pregnancy in combination with parental immunological activation, α-diversity indices based on species richness (number of observed OTUs), estimated species richness (Chao Index) and species diversity (Shannon Index and inverse Simpson) were analysed in linear mixed effect models 37 .  Species richness (numbers of OTUs) was significantly higher on gonadal eggs (GO) compared to brood pouch tissue of non-pregnant males (NP) and on pipefish embryos of early (EP) and late pregnant (LP) males (LMER-#OTU: F 4,67 = 3.84, p = 0.007**; TukeyHSD: GO > NP; GO > EP; GO > LP, Table 1, Fig. 3A). Estimated species richness index (Chao-index) significantly differed among developmental stages (LMER-Chao-dev. stage: F 4,67 = 3.93, p = 0.006**, Table 1). Both gonads and late pregnancy pouch samples revealed a high estimated species richness index (Fig. 3B), whereas pouch tissue of early pregnant and brood pouch tissue of non-pregnant males were in a graduate decline (TukeyHSD-Chao-developmental stage: GO > NP; GO > EP; EP < LP, Table 1, Fig. 3B). We detected an interaction between developmental stage and parental bacteria treatment (LMER-Chao-dev. stage x bact. treat: F 4,67 = 2.53, p = 0.04*, Table 1). Gonad samples of injected females (I) and pouch tissue of injected pipefish males in late pregnancy stage both revealed a significantly higher estimated species richness in comparison to control animals without bacteria injection (C) (TukeyHSD-Chao-dev. stage x bact. treat: GO:C < GO:I; LP:C < LP:I, Table 1, Fig. 3E). In addition, an influence of developmental stage was detected on species diversity indices (LMER-InSimpson: F 4,67 = 4.23, p = 0.004**; LMER-Shannon: F 4,67 = 5.1, p = 0.001**, Table 1, Fig. 3C,D). Bacteria diversity (Shannon index) in pouch tissue of late pregnant males was lower than in pouch tissue of non-or mid pregnant males and in the gonads (TukeyHSD-Shannon: GO > LP; MP > LP; NP > LP, Table 1, Fig. 3C).
Bacterial β-diversity within and between developmental stages and the influence of parental immunological activation. To examine differences in microbiota, we evaluated β-diversity with a Bray-Curtis matrix based on 1419 OTUs obtained by tools implemented in MOTHUR (subsamples > 6000). Bray-Curtis matrix was established on the ratio of shared and unique species relative to the total number of species abundance. We measured variance explained by 'developmental stage' and parental 'bacterial treatment' by applying a permutational multivariate analysis of variance (perMANOVA) in which 'tank' was included as random term by stratifying permutations 10000 times. The bacterial community structure was significantly different between specific developmental stages of gonads and brood pouches (perMANOVA-Bray-Curtis: stage: F 4,67 = 6.26, p < 0.001***, R 2 = 0.24, Table 2). Overall, developmental stage explained 15-24% of total variance assigned to community structure and composition. In contrast, the parental bacteria treatment as single factor did not significantly influence microbial community composition (perMANOVA-Bray-Curtis: bac. treat: F 4,67 = 2.19, p = 0.207, R 2 = 0.02, Table 2). Yet, the significant interaction between developmental stage and parental immunological activation ('bacteria treatment') reveals that the parental immune system influenced the gonad and pouch-specific microbiota (PERMANOVA-Bray-Curtis: dev. stage x bac. treat: F 4,67 = 2.03, p = 0.001**, R 2 = 0.08, Table 2).
An ordination of principal coordinate analysis (PCoA) was applied on OTU abundance (Bray Curtis) measurements for the visualization of each developmental stage (Fig. 4). Subsequently, a linear mixed-effect model (type III sum of squares) was calculated for the first two extracted principle coordinates and significant axes were followed by Tukey HSD post-hoc tests to attain pairewise comparisons between specific treatment groups. In the two-dimensional PCoA depiction (Fig. 4A), the gonads revealed significantly different bacterial communities, they are the most distinct cluster that is set apart along the first principle component (explains 15% of variance) and cluster opposite from all other groups (LMER-Bray-Curtis Axis1: F 4,67 = 17.1, p < 0.001***; TukeyHSD: GO vs NP, EP, MP & LP, Table 3, Fig. 4A,B). Likewise, late pregnancy brood pouch microbiota significantly sets apart along the second principle component (explains 10% of variance), yet it revealed an enormous variance by forming a large ellipse around its centroid rather than a homogeneous cluster (LMER-Bray-Curtis Axis2: F 4,67 = 17.1, p < 0.001***; TukeyHSD: LP vs GO, NP, EP & MP, Table 3, Fig. 4). In contrast, the microbiota of brood pouch  Table 1. Linear mixed effect model of α-diversity indices. LMER was performed on α-diversity indices (average OTUs, Chao, Inverse Simpson, Shannon) including the fixed factors parental bacteria treatment (Bac.treat), developmental stage (Dev. stage), their interaction effect and the random term (tank). Significant LMERs (P < 0.05, indicated in bold letters) were followed by lsmeans (glht function of multcomp package) as post hoc test that includes tank structure. Levels of parental bacteria treatment were: treatment control (C) and injected with heat-killed bacteria (I); Levels of developmental stages were: Gonads (GO), Non-pregnant (NP), early pregnancy stage (EP), mid (MP), and late pregnancy stage (LP) their interaction effect is indicated with a ":" symbol.
SCiEntiFiC RePORTS | (2019) 9:3 | DOI:10.1038/s41598-018-37026-3  tissue of non-pregnant males (NP), early (EP) and mid (MP) pregnant males revealed overlapping centers of gravity and were not significantly different from each other ( Table 3, Fig. 4). We found a developmental stage x parental bacteria treatment interaction effect on the second axis of the PCoA that displays the OTU abundance (LMER-Bray-Curtis Axis2: F 4,67 = 12.6, p < 0.001*** Table 3, Fig. 4). The PCoAs indicate a different ß-diversity structure in the brood pouch during late pregnancy when males had received an immunological activation (Fig. 4). Only late pregnancy brood pouch samples of animals exposed to parental bacteria treatment (Late pregnancy: Injected 'LP:I') showed a significantly different microbiota in comparison to the respective control group (Late pregnancy: Control 'LP:C') ( Fig. 4). To this end, the microbiota on the embryo in late pregnant immunologically activated males (LP:I) formed the most dissimilar cluster that is set apart along the second principle component (LMER-Bray-Curtis Axis2: TukeyHSD: LP:I vs LP:C and all other groups, Table 3, Fig. 4).   Microbiota community compositions according to developmental stage and parental immunological activation. Indicator species analyses 38 and biplot were performed to identify and visualize associations between most abundant microbial species (OTUs) and particular pipefish developmental stages. The biplot consist of a principle component analysis that visualizes a differential clustering pattern with a superimposed factormap demonstrating the contribution of variance (importance %) retained by each bacteria species (Fig. 5).
Results of indicator species analysis conducted using all OUT's can be found as Supplemental Table S3, however, we limited the interpretation of data on the most abundant 50 OTUs. Among those we identified five significant OTUs as indicators for gonadal eggs and three for the late pregnancy stage (Table 4). On gonadal egg surface Brevinema bacteria were overrepresented (Brevinema Otu00044, Otu00031, Otu00014), but also Lewinella (Otu00038) and unclassified_Saprospiraceae (Otu00020) were highly abundant indicators species (Table 4, Fig. 5).
We detected 19 specific OTUs as indicators of the brood pouch microbiota, i.e., the tissue of non-pregnant males and surface of eggs and embryos inside the paternal brood pouch during early, mid and late pregnancy ( Table 4, Fig. 5). They belong to Proteobacteria such as Rhodospirillaceae_unclassified (Otu00007), Acinetobacter (Otu00035), Ruegeria (Otu00004), Marinomonas (Otu00008), Porticoccus (Otu00030), and Phaeobacter (Otu00013) (see Table 4 for further species, Fig. 5). These might shape the microbiota inside the paternal brood pouch that incorporates the embryos during male pregnancy.
We performed a species indicator analysis for developmental stage in combination with parental bacteria treatment (interaction effect) and identified the following indicator species for gonads of females that received an immunological activation: Brevinema (Otu00044, Otu00031), Lewinella (Otu00038) and unclassified_Saprospiraceae (Otu00020) and unclassified_Flavobacteriace (Otu00033) and Sulfitobacter (Otu00026) ( Table 4). The latter two species only occurred in combination with parental bacteria treatment. Also for late pregnancy stage of bacteria treated males, we found overlapping indicator species Kiloniella (Otu00018), Aquimarina (Otu00016), and Ulvibacter (Otu00024) but as a new species Pseudomonas (Otu00034, Otu00042) ( Table 4).

Discussion
In order to understand how microbes colonize reproductively important tissues, persist, change and may outcompete each other to finally shape offspring development and life history, we applied a 16S rRNA gene sequencing approach to identify and compare routes of initial microbial colonisation of embryos. We characterized the parental microbial community over the development of the egg in the gonads and embryo in the brood pouch of the pipefish Syngnathus typhle and assessed the influence of the immune system on microbial colonization during pregnancy by parental bacteria injections. The unique male pregnancy permits to assess maternal contribution  through the eggs and paternal contribution during embryonic development separately. This sheds light on the interaction of microbial colonization and the parental immune system. The most abundant bacteria in the maternal gonads and in the paternal brood pouch of S. typhle belong to Marinomonas (28.0%), Halomonas (10.3%), Aeribacillus (6.3%), Ruegeria (5.8%), Bacteroidetes (4.5%) and Nesterenkonia (3.8%). They are all known members of fish microbiota [39][40][41][42][43] . Marinomonas are important initial colonizers 41,44,45 , Halomonas are gastroinstestinal bacteria 43 and Ruegeria are prevalent on cod larvae 39 .
During mammalian pregnancy the vaginal bacterial community shifts naturally in its structure with respect to diversity and richness 46,47 . The major shift in microbiota during late pregnancy in the pipefish S. typhle suggests the parallel evolution of a similar pattern in male pregnancy. The human vaginal microbiome of pregnant women was suggested to be distinct from non-pregnant women in terms of a decreased species diversity and due to an absence of occasionally present unique taxa 47 . The development of the pipefish brood pouch microbiota community goes in line with this observation: in the brood pouch of the pipefish S. typhle bacterial α -diversity decreased in late pregnancy, while the number of bacteria OTU identified (species richness) was higher. This enhanced variation in the composition of the microbiota suggests a major restructuring that probably facilitates the initial colonisation of the embryos. The passage of surface bacteria into the gut and early digestive tract has been shown to start in oviparous fish when larvae begin to ingest liquid medium 48,49 . Freshly hatched eggs in the brood pouch of the male pipefish (early pregnancy and mid pregnancy) are supplied with essential ions, oxygen and nutrition proteins over a placenta-like structure from the father in addition to the essential supply by the maternal yolk sac [50][51][52][53] . The embryonal development takes place in a paternally shaped environment. To this end, the embryos are supposedly soaked in a cocktail of paternal bacteria. This makes it tempting to speculate that once yolk sac supply is depleted, which usually corresponds to the opening of the larval mouth (i.e., between mid and late pregnancy in pipefish), the paternal microbial community are dominating the gut of pipefish larvae.
In mammalian viviparity the embryo is supposed to be kept almost sterile in the mother's womb, in pipefish the environment could play a pronounced role in the microbial colonisation of the developing embryos. Due to  88 . Asterix symbol behind OTU number * , ** , *** demonstrate significant Indicator species according to indicator species analysis (  excessive larval growth towards the end of pregnancy and the forthcoming birth, the skin fold that closes the brood pouch upon intake of the eggs, gets chapped, which increases the permeability of the brood pouch. The embryos thus get in contact with the bacteria of the natural habitat (surrounding water). This invasion of environmental microbes could prime or pre-adapt the embryos towards the microbial community, they will encounter in the near future. The shift from the constant microbiota during early and mid-pregnancy in the pouch towards a more variable microbiota during late pregnancy could support the hypothesized invasion of environmental microbes. Drivers of the microbial modulation are Kiloniella, Aquimarina and Ulvibacter. Kiloniella have been previously reported to increase in abundance on the skin of fishes in response to stress 40 , while Aquimarina not only belong to fish gut microbiota but also have the potential to inhibit the growth of pathogenic bacteria 62 . Ulvibacter belong to Flavobacteria that are common bacterial colonizers in the gastrointestinal tract of marine fishes 63 . A shift in the bacterial community during late pregnancy is further substantiated by a common pool of indicator species shared over all developmental stages, with the exception of the late pregnancy stage. Due to high abundance (under the top 20 most abundant OTUs), these bacteria might be important for the early embryonic development. Among them are Caldalibacillus, Virgibacillus and Pseudoidiomarinae, Nesteronkonia, Bacillaceae and Phyllobacterium. Virgibacillus has been reported to have antimicrobial activity, which plays an important role in host defence against pathogenic bacteria in the ovaries or in the pouch 64,65 . These microbes potentially shape initial colonization that occurs prior to the mouth opening, i.e. transovarially, or they are relocated with the paternal nutrient transfer in the pregnancy to the embryos.
This study provides first insight into potential maternal and paternal contribution to offspring microbiota. Syngnathids are an enigmatic model system to assess this, as only the evolution of male pregnancy permits a straightforward differentiation between microbial colonization of the ovary from the colonization during pregnancy. As females produce the eggs and males brood the embryos in their brood pouch, our data simultaneously reflect sex-specific microbial contribution. We identified a microbial community on the eggs in the maternal ovaries that was distinct from all other developmental stages. The bacteria driving this difference are associated to fish mucus (Brevinema) 66 and fish gut microbiota (Brevinema, Saprospiraceae) 67,68 .
While highly prevalent in the gonads and on the pipefish eggs, none of these bacteria was abundant during late pregnancy. Promising candidates for a transovarial transfer are thus Sulfitobacter and members of the Flavobacteriaceae, as they were prevalent both in the gonads and in the pouch during late pregnancy but not in the pouches of non-pregnant males or during early pregnancy. This makes an involvement in both a transfer via eggs into the brood pouch and in initial microbial colonization likely. Both Sulfitobacter and Flavobacteriaceae are important members of fish larval commensal microbiota 41,58 .
The parental immune system is in close interaction with the microbial community in the gonads and in the parental brood pouch. The identified expansion in species richness during late pregnancy was much more pronounced upon paternal bacterial treatment envoking an immunological activation. This suggests that the immune system of the parents may shape the microbial community in the gonads, the time point when transovarial transfer of microbiota can occur, but also during late pregnancy, when initial microbial gut colonization is most likely to take place. Some of the bacterial indicator species for the gonad and late pregnancy microbial communities are also overrepresented in the microbial community of the gonads from females that were previously injected with heat-killed bacteria (interaction of developmental stage and parental bacteria treatment). In addition to the previously discussed indicator species (Kiloniella, Aquimarina, Ulvibacter) Pseudomonas is highly prevalent during late pregnancy, if the immunological activation is considered. Pseudomonas belongs to the dominant commensal microbiota of marine fish species 58 . Marinomas was most abundant during late pregnancy on the embryos of fathers that were previously exposed to bacteria treatment.
The previously described bi-parental immune priming in the pipefish S. typhle [30][31][32] , which entails that the immune status of both parents shapes the immunological performance of their offspring, could directly interact with the microbial community in the ovaries and in the brood pouch. As such, immunological activation imposed over parental bacteria treatment may change the bacterial composition to higher abundance of Kiloniella, Aquimarina, Ulvibacter & Marinomonas bacteria. These are commensal bacteria with the potential to specifically aid in fighting virulent bacterial infections or alternatively, to boost the offspring immune response.

Conclusion
Our data suggest a co-adaptation between parental investment, the host microbiota and the immune system during pipefish male pregnancy. Immune defences and the microbial community may be simultaneously transferred from mothers and fathers to the offspring, which can have substantial sex-specific developmental implications that need to be assessed in more detail in the future. Such co-adaptation between microbiota and the vertebrate adaptive immune system are likely to not only enable the clearance of potentially virulent pathogens but also shape the persistence of specific mutualistic microbial communities 33 .

Material and Methods
Experimental setup and sampling. Pipefish (females and pregnant/non-pregnant males) were collected from Kiel Fjord (54°44′N; 9°53′E, July 2014) and hosted under Baltic summer conditions (15 PSU;18 °C) in 500 L tanks. Pipefish were separated according to their sex, and males also according to their developmental stage: 1.group: Females; 2.group: Non-Pregnant males (NP); 3.group: Early Pregnant males (EP) (1-2 weeks); 4.group: Mid Pregnant males (MP) (2-4 weeks); 5.group: Late Pregnant males (LP) (4-6 weeks) (see Fig. 1). 20 females, 20 males per pregnancy stage (20 early, 20 mid 20 late pregnancy), and 40 non-pregnant males were used for the study (in total 120 fish). In the following, 10 females and 10 pregnant males/per pregnancy stage (EP/ MP/ LP) as well as 20 non-pregnant males (undeveloped brood pouch) were randomly chosen and injected with 1:1 mixture of two heat-killed bacteria (Vibrio spp. Italy2K3 and Tenacibaculum maritimum) according to 32 to boost their immune system and trigger a parental immune transfer. The bacteria treated individuals and the non-treated pipefish were randomly transferred into 200 L aquaria (3 tank replicates per group (injected vs. non-injected), 20 pipefish per tank). Pipefish were fed twice a day with frozen mysids. All tanks were connected through the same filtered local Baltic Seawater (semi-flow through circulation). After 10 days of incubation, the previously immune-challenged pipefish received a secondary injection with the same dosage of heat-killed Vibrio bacteria. The ten days between the two paretntal bacteria treatements were chosen to permit sufficient time for an impact of the immune system on the microbiota diversity and simultaneously allow for parental transfer of immunity. 20 hours after the second immune challenge, pipefish were killed by an overdose of an anesthetic (MS-222). The gonads of females were dissected and microbial samples were taken with a sterile swab at the surface of the gonadal eggs. The brood pouch of pregnant pipefish males was opened and microbial samples were taken from the surface of the embryos inside the paternal brood pouch with a sterile swab. The brood pouch tissue of non-pregnant males was opened and microbiota samples were taken with sterile swabs. DNA extraction, Library preparation and sequencing. Microbial DNA was extracted with the Dneasy 96 Blood & Tissue Kit (Qiagen) following the manufacturers protocol, with ameliorations according to Franke et al. (2017) 69 . The dual-index sequencing strategy developed by 70 was applied. 16S rRNA genes of samples and a negative (sterile swab) control were amplified (Phusion High-Fidelity DNA Polymerase, Thermo Fisher Scientific) with the labelled primer pairs F515: (AATGATACGGCGACACCGATCTACAC <i5> TATGGTAA TTGTGTGCCAGCMGCCGCGGTAA); R806: (CAAGCAGAAGACGGCAACGAGAT <i7> AGTCAGTCAG CCGGACTACHVGGGTWTCTAAT) spanning the hypervariable variable region V4 of the 16S rRNA gene 71 . Forward and reverse primer pairs contained adapters, barcodes, pad and linker sequences as described by 70 . Several negative controls without sample were included per plate. The PCR was carried out under the following conditions: 98 °C for 30 s, 30 cycles of 98 °C for 9 s, followed by 55 °C for 15 s and 72 °C for 20 s with a final elongation step at 72 °C for 10 min. PCR products were purified with the MinElute 96 UF PCR purification kit (Qiagen). DNA concentrations were measured by spectrophotometry (NanoDrop ND-1000, Peqlab). PCR products were pooled to one subpool per plate in an equimolar concentration (~30 ng per sample), respectively, run on a 2% agarose gel, and amplicons of the expected size (~300 bp) extracted using the NucleoSpin gel and PCR clean-up kit (Macherey-Nagel). The extraction products of subpools were fluorometrically quantified (Qubit fluorometer, Invitrogen) and pooled in equimolar concentrations. The amplicon library was sequenced on the Illumina® MiSeq using the MiSeq Reagent Kit v2 500 cycles sequencing chemistry and mixed with Illumina generated phiX control libraries.
Sequence processing and data analysis. Raw sequences were de-multiplexed using Casava v.1.8.2, assembled and filtered with the software MOTHUR v.1.16.1 according to the MiSeq SOP pipeline 72,73 . Adaptor, tag, and primer sequences were removed from raw sequences. Unique sequence reads were merged and aligned against the SILVA alignment database (release #119) 74,75 . All misaligned sequences not covering the variable region 4 were removed (SILVA alignment position 1968 to 11550) 74,75 . To reduce sequencing noise, a pre-clustering step (2 bp difference) was performed 72 and chimeric sequences were eliminated by the UCHIME algorithm implemented in MOTHUR 76 . The taxonomy of all sequences was estimated using the classify.seqs function in MOTHUR against the Ribosomal Database Project (RDP) training set v.9 77 with a 80% bootstrap threshold. Non-target sequences (e.g. chloroplasts, mitochondria, eukaryotic 18S rRNA) were removed. The sequences were clustered at the 0.03 difference level to obtain operational taxonomic units (OTUs) with average neighbour algorithm. By randomly taking the same number of sequences from each sample we rarefied the samples to 6000 reads each, which enables to perform and compare diversity measurements at the same sequencing depth (Supplemental Table 1 Table S1), specific statistical tests (type III sum of squares, permutational tests) were applied for the analysis. Species richness (number of observed OTUs), estimated species richness (Chao 1 Index, corrected for sample size) and species diversity (Shannon Index and inverse Simpson) were calculated with the summary.single command implemented in MOTHUR based on a dataset subsampled to a number of 6000 reads per sample (Supplemental Table S1). For visualization and interpretation of the microbial community data, we used standardized 97%-OTUs for relative abundance (Bray-Curtis distance matrix) analyses.
Statistical Analysis. All statistical tests and visualizations were performed using a rarefied subset of 6000 reads/sample in the R environment (R Core Team 2015). The α-diversity measurements species richness (number of observed OTUs), estimated species richness (Chao Index), and species diversity (Shannon Index and inverse Simpson) were analysed fitting for each index a linear mixed effect model 34 . For the statistical model we applied the fixed interaction term 'developmental stage' x 'bacteria treatment' , 'tank' was included as random term. LMER models were performed with the lmer function implemented in the lme4 package of R 78 by using type III sum of squares and Satterthwaite approximation for the degrees of freedom. For multiple comparisons of fixed and interaction terms all significant LMERs were followed by post-hoc t-tests with the ghlt function associated in the multcomp package of R 79 .
β-diversity measurements were assessed to analyse bacterial communities between tissues and treatment groups based on abundance (Bray-Curtis) of shared 97%-OTUs with the vegan package v. 2.3-0 in R 80 . We evaluated changes of microbiota at the surface of the eggs in the ovaries, and at surfaces of increasing developmental embryonal stages in combination with parental bacteria treatment by applying a permutational multivariate analyses of variance (perMANOVA) with the adonis function of the vegan package in R 81 . PERMANOVA models were based on abundance (Bray-Curtis) data of all OTUs identified in the samples, applying 'developmental stage' × 'bacteria treatment' as fixed factors and stratifying 10000 permutations within each tank replicate. Visualization of variations in microbiota among developmental stages (gonads, brood pouch of non-pregnant, early, mid and late pregnant males) and parental bacteria treatments were assessed with analysis of principal coordinates (PCoA) based on bray Curtis distance matrix using the pcoa function from the ape R-package [82][83][84] . A PCoA is based on eigenvalue equation but can use any dissimilarity index and distances between points in the plot reflecting original distances 85 . Hypothesis-based treatments were added as dispersion ellipses to the ordination plots with 0.95 confidence intervals. We extracted the principle coordinates (scores) of the first two axes and fitted a linear mixed-effect model for each single axis by applying the lmer function implemented in the lme4 package of R 86 . The scores of the first two principle components were used for the statistical analysis to attain the projection that accounts for the most relevant variation 87 . We fitted a linear mixed-effect model using 'developmental stage' and 'bacteria treatment' as fixed factors and 'tank' as random term. The linear mixed-effect model was performed with type III sum of squares and Satterthwaite approximation for degrees of freedom 86 . Significant axes were followed by Tukey HSD post-hoc test with the ghlt function associated in the multcomp package of R 79 to attain pairewise comparisons between specific treatment groups. Bacterial distribution patterns and diversity were based on the 50 most abundant OTUs that contribute 90% of all sequence counts (Indicator species on all OTUs can be found in Supplemental Table S3). Indicator bacteria species and associations between species (OTUs) and tissue of interest (gonads, brood pouch of non-pregnant, early, mid and late pregnant males in combination with parental bacteria treatment vs. no bacteria treatment), were identified with the package indicspecies in R based on 10000 permutations 38 . By drawing a biplot (factor map) with the factoextra package implemented in R 88 in which a colour gradient highlights most important species (OTUs), we visualized the contribution of variance (%). In addition, similarity-clustered heatmaps of the 50 most abundant OTU species were conducted. Heatmaps were illustrated with the aheatmap function of the NMF package in R by applying an "euclidean" distance and complete linkage method for hierarchical clustering.
Ethics approval and consent to participate. All animals were handled according to the animal welfare laws of Germany, under a permit of the "Ministerium für Landwirtschaft, Umwelt und ländliche Räume des Landes Schleswig Holstein" called "Komparative Vergleichsstudie von Immunantworts-Transfer von Eltern zu Nachkommen in Fischarten mit extremer Brutpflege".