In this study, epiphytic and surrounding seawater bacterial samples were collected at six developmental time points from mature sporophytes, sporelings to juvenile sporophytes stages at Qingyutan, Rongcheng, Shandong Province, China (37°10' 2.19" N, 122°34' 54.8" E). Mature sporophytes (MS) with sporangia were collected on 8^th August, 2019 and were rinsed with sterile seawater to remove the loose epiphytes on the surface ( Figure 1 ). Epibacteria were sampled by swabbing 50 cm² of the sporangia surface with sterile cotton swabs near the one-third of mature sporophytes away from the meristematic tissue. The swabs with epibacteria were stored in 50 mL sterile Eppendorf tubes. Usually it takes one month for the zoospores to develop into sporelings (S), which cover the substrate curtains (40 ×115 cm) in the greenhouses. The environmental factors, such as temperature (10 ± 1 °C), salinity (30‰) and light intensity (80 μE m^-2 s^-1), were controlled in the greenhouses. Method for sampling epibacterial samples of sporelings referred to Han et al. (2021). Epibacteria of 4-week-old sporelings (S1, 0.1-0.2 mm in length), 7-week-old sporelings (S2, 1-5 mm in length) and 9-week-old sporelings (S3, 10-20 mm in length) were sampled by swabbing on the substrate curtains by using sterile cotton swabs on 7^th September, 27^th September and 11^th October 2019, respectively ( Figure 1 ). Swabs with epibacteria were put into 50 mL sterile Eppendorf tubes. Juvenile sporophytes (JS) are referred to those sporelings which were transferred into the sea and cultivated for a period of 1-2 month. Epibacterial samples for both 4-week-old juvenile sporophytes (JS1, 30-40 cm in length) and 9-week-old juvenile sporophytes (JS2, 50-60 cm in length) were collected on 29^th November 2019 and 4^th January 2020, respectively ( Figure 1 ). The sampling method of epibacteria on juvenile sporophytes were the same as mature sporophytes.

Regarding to the collection of seawater bacterial samples, 2.5 L surrounding seawater were collected by using sterile bottles at each developmental time points. Seawater samples were firstly filtered through the membrane of 3 μm pore size to wash off the larger biotic or abiotic particles, and then were filtered through a membrane of 0.22 μm pore size (Whatman, UK). These membranes with seawater bacteria were put into 50 mL sterile Eppendorf tubes. All bacterial samples were kept frozen at -80 °C until DNA was extracted. Six replicates for both epibacterial and seawater bacterial samples were designed for each developmental time points, respectively.

We used HiPure Soil DNA kits (Magen, Guangzhou, China) to extract DNA of bacteria following the manufacturer’s instructions. The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified by PCR (95°C for 2 min, followed by 27 cycles at 98°C for 10 s, 62°C for 30 s, 68°C for 30 s and a final extension at 68°C for 10 min) with the primers 341F (5′- CCTACGGGNGGCWGCAG-3′) and 806R (5′- GGACTACHVGGGTATCTAAT -3′). PCR reactions were performed in 50 μL mixture, containing 10 × Buffer KOD (5 μL), dNTPs (5 μL, 2 mM), primers 341F and 806R (1.5 μL of each, 10 μM), KOD DNA Polymerase (Toyobo, Japan) and 100 ng template DNA.

PCR amplicons were purified from 2% agarose gels using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). Quantification was performed using ABI StepOnePlus Real-Time PCR System (Life Technologies, Foster City, USA). The 2 × 250 bp paired-end reads were generated with Guangzhou Genedenovo Biotechnology (China) on the Illumina Novaseq6000 platform (Guo et al., 2017).

The DADA2 (version 1.14) was used to denoise and remove low-quality reads (Callahan et al., 2016). Then paired end denoised reads were merged as tags with a minimum overlap of 12 bp. Chimera sequences were identified and were deleted by UCHIME algorithm (Edgar et al., 2011). The high-quality tags were clustered into ASVs (amplicon sequence variants) with 100% similarity. The confidence threshold values ranged from 0.8 to 1.0. All ASV sequences were classified into organisms by a naive Bayesian model using RDP classifier 2.2 (Wang, 2007) based on SILVA database (version 132) (Pruesse et al., 2007). Because there were much higher standard errors in some groups, we removed the three replicates with outliers. Thus, finally, only three replicates close to the average value were selected for the analysis in this study. The sequences of chloroplast, mitochondria and archaea origin were removed from the dataset.

Genera present in all replicates of epiphytic or surrounding seawater bacterial communities were selected and defined as the core genera. Core genus analysis of both epiphytic and surrounding seawater bacterial communities were performed in R project Venn Diagram package 1.6.16 (Chen and Boutros, 2011). Further, top core genera were identified by screening the core genera with relative abundance more than 1.0%, and was visualized in R project ggplot2 package 2.2.1 (Wickham, 2011).

The stacked bar plot of the communities composition was visualized in R project ggplot2 package 2.2.1 (Wickham, 2011). The Chao1, Pielou’s and Shannon diversity indices were calculated in QIIME (version 1.9.1) and the comparison between samples were calculated by Welch’s test. Non-metric multidimensional scaling (NMDS) and Permutational multivariate analysis of variance (PERMANOVA) were calculated based on unweighted unifrac dissimilarities in R project Vegan package 2.5.3.

Indicator species were those genera which influence structure differences in both epiphytic or surrounding seawater microbial communities. Indicator species were analyzed by calculating the indicator value (IndVal) in labdsv package (version 2.0-1) of R project and performing cross-validation test (Roberts, 2016) based on the abundance and occurrence frequency at genus level. Genera with IndVal ≥ 0.7 and p value ≤ 0.05 was selected as indicator genus (Glasl et al., 2019).

There were 18 epibacterial and 18 seawater bacterial samples were analyzed from mature sporophytes, sporelings to juvenile sporophytes of northern farmed S. japonica ( Figure 1 ). A total of 4,653,122 paired-end raw reads and 4,499,030 clean reads were obtained after filtering low quality reads and removing the mitochondria, chloroplasts and eukaryotes sequences. The 24,821 amplicon sequence variants (ASVs) and were assigned across 36 samples. Among the samples, 14,035 ASVs were clustered in the epibacterial (EB) group, which much higher than those in the seawater (SW) group (13,406 ASVs). Only 10.6% ASVs were shared between epiphytic and seawater bacterial communities ( Supplementary Figure 1 ). ASVs were classified into 41 phyla, 111 classes, 247 orders, 372 families and 958 genera. The rarefaction curves showed good diversity coverage (> 99.2%), indicating the sequencing amount was sufficient ( Supplementary Figure 2 ).

Temporal shifts of epibacterial communities associated with northern farmed S. japonica were analyzed from mature sporophytes to juvenile sporophytes ( Figure 2 ). In the epibacterial communities, Proteobacteria (mean relative abundance: 49.1% to 71.1%) was the most dominant phylum at each time point ( Figure 2A ). At the genus level, Alcanivorax was the most dominant genus at both MS-EB (20.4%) and S1-EB (5.5%). Bacillus was the most dominant at S2-EB (19.9%) and S3-EB (17.0%). Halomonas and Cobetia were the most dominant genera at JS1-EB (12.1%) and JS2-EB (15.5%), respectively ( Figure 2B ). Regarding to the surrounding seawater bacterial communities, the most dominant phylum was Firmicutes at MS-SW (85.8%) and S3-SW (33.5%). Actinobacteria was the most dominant phylum at S1-SW (51.5%). Meanwhile, Proteobacteria was the most dominant phylum at S2-SW(42.2%), JS1-SW(71.8%) and JS2-SW (52.5%). From mature sporophytes, sporelings to juvenile sporophytes, the relative abundance of Firmicutes decreased from 85.8% to 24.6%, and Proteobacteria increased from 6.5% to 62.1% ( Figure 2A ). At genus level, the most dominant genera were Virgibacillus at MS-SW (30.8%), Rubritalea at S2-SW (29.2%), Stenotrophomonas at JS1-SW (22.0%) and Cobetia at JS2-SW (25.2%). Oceanobacillus was the most dominant genus at both S1-SW (4.4%) and S3-SW (19.6%) ( Figure 2B ).

There were significant differences in Chao 1 indexes among six developmental time points (p < 0.05). Except S1-EB, the Chao 1 (species richness estimates) indexes of epibacterial communities exhibited increasing trend with the development with northern farmed S. japonica ( Figure 3A ). Regarding to Pielou’s evenness and Shannon indexes, there were no significant shifts with the development of S. japonica ( Figures 3B, C ; p > 0.05). In surrounding seawater bacterial communities, Chao 1 indexes were significantly different with development of S. japonica ( Figure 3D , p < 0.05). The Pielou’s evenness was not significantly different ( Figure 3E , p > 0.05), but the Shannon indexes exhibited significant differences ( Figure 3F , p < 0.05). At the same developmental time point, the Chao 1 indexes ( Table 1 , p < 0.01) and Pielou’s evenness ( Table 1 , p < 0.05) of epiphytic and seawater bacterial communities were significantly different at MS and JS2. While, the Shannon indexes of epiphytic and surrounding seawater bacterial communities exhibited significant differences at JS2 ( Table 1 , p < 0.05).

The beta diversity of epibacterial communities displayed significant differences at six developmental time points ( Figure 4A ). PERMANOVA showed that developmental time explained 40.8% of epibacterial communities variation (PERMANOVA: F value = 1.66, R^{2 =} 0.41, P = 0.001, Table 2 ). Meanwhile, there were significant differences in the beta diversity of seawater bacterial communities at the different time points (PERMANOVA: F value = 1.89, R^{2 =} 0.44, P = 0.001, Table 2 ; Figure 4B ). Moreover, there was no significant difference between the structure of the epiphytic and the surrounding seawater bacterial communities. (PERMANOVA: F value = 1.2, R^{2 =} 0.03, P = 0.058, Table 2 ; Supplementary Figure 3 ). The structure of epiphytic and surrounding seawater bacterial communities were also not significantly different at the same developmental time point by pairwise comparisons ( Supplementary Table 1 , p > 0.05).

In order to further explore the structure differences of epiphytic or surrounding seawater bacterial communities at six developmental time points, indicator species were analyzed based on IndVal ≥ 0.7 and p value ≤ 0.05 at genus level. The indicator species of MS-EB were Gracilimonas, Vibrio, C1-B045. Thalassospira was indicator species at the S1-EB. Whereas, Pseudomonas Rhodopirellula and Salegentibacter were the indicator species at JS1-EB. There were no indicator species at S2-EB, S3-EB and JS2-EB ( Figure 4C ). In addition, bubble plot showed different indicator species at different developmental time point to reflect the structure of temporal differences in surrounding seawater bacterial communities ( Figure 4D ). Virgibacillus and Sediminibacillus were indicator species at the MS-SW. Rubritalea was indicator species at the S2-SW. Acinetobacter, Brevundimonas and Pseudomonas were indicator species at the JS1-SW.

In this study, core genera were referred to those that were present in all samples at the developmental time points, regardless of their relative abundance. Overall, there were 48 core genera in the epibacterial communities ( Supplementary Figure 4A ). Among these core genera, the mean relative abundance greater than 1.0% were referred to top core genera. The top core genera included Bacillus, Alcanivorax, Halomonas, Psychrobacter, Blastopirellula, Erythrobacter, Rubritalea, Pseudomonas, SM1A02, Muricauda and Gimesia ( Table 3 ). From MS to S3, the relative abundance of Bacillus in the epibacterial communities gradually increased (ranging from 0.1% to 17.0%) ( Supplementary Figure 5 ). There were 29 core genera ( Supplementary Figure 4B ) and 7 top core genera ( Table 3 ) in the seawater bacterial communities. Among these top core genera, the relative abundance of Bacillus (27.2%), Virgibacillus (31.9%) and Sediminibacillus (15.8%) were higher at MS-SW. While, the relative abundance of Bradyrhizobium (4.6%) and Rubritalea (23.4%) was higher at both S2-SW ( Supplementary Figure 5 ). Acinetobacter (10.4%) and Pseudomonas (10.4%) was higher at both JS1-SW. Bacillus, Rubritalea and Pseudomonas were top core genera between epiphytic and surrounding seawater bacterial communities. From MS to S2, the relative abundance of Bacillus in the surrounding seawater bacterial communities gradually decreased, which was the opposite of the epibacterial communities ( Supplementary Figure 5 ). The top core genera Rubritalea between epiphytic and surrounding seawater bacterial communities had the highest relative abundance at S2, and the highest relative abundance of Pseudomonas at JS1. The relative abundance of Rubritalea and Pseudomonas were higher in the surrounding seawater bacterial communities at S2 and JS1 ( Table 3 ).

It was widely accepted that dominant taxa in microbial communities of seaweeds were species- specific at genus level (Bengtsson and Øvreås, 2010; Mancuso et al., 2016; Lemay et al., 2018b, Lemay et al., 2021a; Guo et al., 2022). In this study, the dominant genera with the highest relative abundance were Alcanivorax (MS-EB and S1-EB), Bacillus (S2-EB and S3-EB), Halomonas (JS1-EB) and Cobetia (JS2-EB) ( Figure 2B ). Previous studies had shown that Flavobacterium, Yoonia- Loktanella and Planctomyces were the most dominant genera for brown alga Sargassum horneri (Mei et al., 2019), Sargassum thunbergii (Guo et al., 2022) and the farmed Caulerpa lentillifera (Pang et al., 2022), respectively. This difference indicated that dominant genera were seaweeds specific (Roth-Schulze et al., 2016; Korlević et al., 2021). What’s interesting, in this study, the dominant genera Alcanivorax, Bacillus and Halomonas were also the core genera. Alcanivorax is a degrading bacterium for marine hydrocarbon pollutants (Olivera et al., 2009), which may render its seaweed hosts to decontaminate the hydrocarbon pollutants and contribute to the healthy development. Bacillus, especially Bacillus subtilis, has strong antibacterial activity against the pathogens Aeromonas hydrophila, Vibrio vulnificus and Vibrio parahaemolyticus of the brown alga Sargassum myriocystum (Chakraborty et al., 2017). Moreover, Halomonas has been found at the wide range of pH and temperatures as well as at almost any range of salinity, which enables Halomonas to colonize the variety of habitats (Kim et al., 2013). Halomonas is also beneficial to the morphogenesis and development of seaweeds, which excretes a specific regulator like cytokinins that enhanced cell division (Spoerner et al., 2012). For example, Halomonas sp. Z3 has been reported to promote the development and increase the number of individuals of the brown alga Ectocarpus sp. (Tapia et al., 2016). Based on the relative abundance advantage and beneficial effects, we speculate that Alcanivorax, Bacillus and Halomonas may play beneficial roles during the development of northern farmed S. japonica. However, their effects on the seaweed development still remain to be further verified by isolating the pure bacterial strains.

The temporal shifts of diversity with the development of seaweeds was one of the evidences of microbial communities shifts (Weigel and Pfister, 2019). Previous studies indicate that the alpha diversity gradually increased with the development of Laminaria saccharina (Staufenberger et al., 2008), Cystoseira compressa (Mancuso et al., 2016) and the Sargassum muticum (Serebryakova et al., 2018). Consistent with these previous researches, we found the Chao 1 indexes of epibacterial communities increased and was significantly different with the development of northern farmed S. japonica (p < 0.05). Meanwhile, beta diversity revealed significant differences in the epibacterial communities among different developmental stages ( Figure 4A , p < 0.01). The epiphytic bacterial communities exhibited a significant increase in species diversity as well as significant differences in community structure, but the Pielou’s evenness was stable during the development of northern farmed S. japonica ( Figure 3B ). This suggested that the epibacterial communities continued to recruit new bacteria to assemble the communities with the development of S. japonica, rather than based on the increase in the abundance of already colonized bacteria. In addition, it is reported that different metabolites (e.g., enzymes and phenolics) can be secreted by the seaweeds at different development time points (Collén and Davison, 2001; Malik et al., 2020; Han et al., 2021; Lemay et al., 2021a), which can attract the attachment of certain bacteria (Zheng et al., 2005) and thus have the unique microbial composition and structure at different developmental time points.

Studies have shown that epibacterial communities of seaweeds were specific and were influenced by the surrounding seawater bacterial communities (Fahimipour et al., 2017; Lemay et al., 2018a; Cleary and Huang, 2020; Juhmani et al., 2020). It was found that the dominant bacteria of Macrocystis pyrifera (Michelou et al., 2013), Thalassia testudinum and Syringodium filliforme (Ugarelli et al., 2018) were shared between the epiphytic and surrounding seawater communities at the phylum level. Consistent with these previous investigations, we observed that Proteobacteria was the dominant phyla in both epiphytic and surrounding seawater bacterial communities. Moreover, the top core genera in the epibacterial communities were all present in the surrounding seawater bacterial communities, which is similar to the results in Nereocystis luetkeana and M. pyrifera (Weigel and Pfister, 2019). The top core genera Rubritalea and Pseudomonas between epiphytic and surrounding seawater bacterial communities had the highest relative abundance at the same developmental time point, and their relative abundance were higher in the surrounding seawater bacterial communities. Meanwhile, the structure of the epiphytic and surrounding seawater bacterial communities did not differ significantly. These results suggested that the surrounding seawater bacterial communities influenced the epibacterial communities. It could be explained by the fact that the surface of S. japonica was constantly in contact with the diversified free-living bacteria in the surrounding seawater and most of the epibacteria are recruited from the surrounding seawater (Mancuso et al., 2016; Weigel and Pfister, 2019). It is worth noting that microbial communities of seaweeds have been shown to be influenced by a variety of factors, including host species, geography and environment, and so on (Hollants et al., 2013; Morrissey et al., 2019). Therefore, the factors influencing the shifts of epibacterial communities of the farmed S. japonica will need to be explored more fully and deeply in the future.

It has been found that the epibacterial communities associated with seaweeds are also influenced by the seaweeds themselves (Weigel and Pfister, 2019; James et al., 2020; Lemay et al., 2021b). In this study, the alpha diversity of the epibacterial communities, including Chao1, Pielou’s evenness and the Shannon indexes, were significantly higher than those of surrounding seawater bacterial communities at JS2 (p < 0.05). The reason for the higher alpha diversity at JS2 maybe due to the faster growth of S. japoninca. Zimmerman and Kremer (1986) and Zhang et al. (2020a) reported that the concentrations of mannitol and laminarin on S. japonica reduced when it grew fast. Therefore, we speculate that the decreased concentrations of mannitol and laminarin on S. japonica at JS2 may attract or inhibit bacterial attachment and then led to a higher alpha diversity of the epibacterial communities. This indicated that the epibacterial communities associated with S. japonica were also influenced by the seaweeds themselves. Our results were consistent with previous findings obtained from Ecklonia radiata (Marzinelli et al., 2015), Cystoseira compressa (Mancuso et al., 2016) and Caulerpa taxifolia (Morrissey et al., 2019). This suggested that the epibacterial communities were more affected by the host itself than the surrounding seawater bacterial communities at the faster growth stage of S. japonica.