A total of 460 undamaged and uniformly sized fertile eggs of Ross 308 commercial broilers were purchased from the Shuncheng Broiler Breeder Farm (Ningxiang, Hunan). All eggs were collected from 35-week-old flocks. The animal procedures in this study were approved by the Institutional Animal Care and Use Committee of Hunan Agricultural University (GBT2018). The fertile eggs were randomly divided equally into group CS and group CC and the initial weight was recorded. There was no significant difference in the initial egg weight between the two groups of fertile eggs (CS: 68.37 ± 3.75 g, CC: 67.81 ± 4.17 g, p = 0.135, n = 460). Before incubation at 37 ± 0.5°C with 60 ± 5% relative humidity, the shell surface of group CS was sterilized with 70% alcohol, while that of group CC was not sterilized. The alcohol was sprayed evenly over the entire surface of the egg and then wiped clean with sterile cotton. After sterilization of each egg, the gloves were sterilized with alcohol, followed by the surface sterilization of the eggshell of the next egg. A sterile marker (DOTCH Puru Marker) was used to number all eggs. The incubators were sterilized by ultraviolet light sanitization before use, and the two groups were separated. Eggs were transferred to pedigree hatch baskets at E18. Each egg was isolated from the others in individual cells in the hatcher tray, thus retaining the treatment information of each chick upon fertilization. All procedures were performed in a sterile room with all the manipulations in a sterile hood. The ventilation system was connected to the common environmental air supply without any sterilizing treatments. All the eggs were candled at E05 to check the fertility and the unfertilized eggs were eliminated. In total, 24 fertile eggs from each of the two groups were randomly selected to collect samples (E0, E07, E15, and E21) and to weigh the embryos at the indicated time points (E07, E15, and E21) (3), respectively. Chicken embryo growth performance was measured by calculating the ratio of chicken embryo weight to initial egg weight. In total, 12 of these 24 eggs were selected to investigate the microbiota in the egg yolk. The eggshell was gently peeled off using a sterile tweezer to expose the embryo, and then, the yolk samples (5 mL) were collected by puncturing the yolk sac membrane with a syringe and homogenized. Subsequently, the abdomens of chicken embryos were dissected with a scalpel and then the whole gut samples were carefully collected on a 4°C sanitized working bench. The yolk samples and gut samples were stored at –80°C after snap freezing in the liquid nitrogen. After hatching, the number of hatched chickens was recorded to determine the hatchability with the following formula: Hatchability = (No. of hatched eggs/No. of fertile eggs after sampling) × 100.

Approximately 5 mL of egg yolk samples were thawed on ice, then put in 10 mL of sterile 1× PBS buffer to remove the excessive yolk fat, and centrifuged at 16000 rpm for 5 min, and DNA was extracted from yolk samples using the TIANamp Micro DNA Kit (TIANGEN, cat#DP316). Total DNA was extracted from homogenized embryonic intestinal samples using the TIANamp Stool DNA Kit (TIANGEN, cat#DP328), according to the manufacturer’s instructions. The concentration and quality of the isolated DNA were assessed by using a NanoDrop spectrophotometer (Thermo Scientific, United States). Amplicons of the V4 and V5 hypervariable regions of 16S rRNA were amplified by using the sample-specific sequence barcoded fusion primers (forward 5’-GTGCCAGCMGCCGCGGTAA-3’ and reverse 5’-CCGTCAATTCMTTTRAGTTT-3’). Sample-specific 7-bp barcodes were incorporated into the primers for multiplex sequencing. The PCR components contained 5 μL of Q5 reaction buffer (5×), 5 μL of Q5 High-Fidelity GC buffer (5×), 0.25 μL of Q5 High-Fidelity DNA Polymerase (5 U/μL), 2 μL (2.5 mM) of dNTPs, 1 μL (10 uM) of each forward and reverse primer, 2 μL of DNA template, and 8.75 μL of ddH₂O. Thermal cycling consisted of an initial denaturation at 98°C for 2 min, followed by 35 cycles consisting of denaturation at 98°C for 15 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s, with a final extension of 5 min at 72°C. PCR amplicons were purified with Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN) and quantified using the PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). 16S rRNA sequencing was performed using the Illumina Novaseq_PE250 (Illumina) sequencing platform. Sequencing services were provided by Personal Biotechnology Co., Ltd., Shanghai, China. The data were analyzed by using the free online platform Personalbio GenesCloud¹ and a software package for analyzing taxonomic or metabolic profiles (STAMP v2.1.3).

In total, 13,356,685 high-quality reads were generated for further analysis. The Quantitative Insights Into Microbial Ecology (QIIME2, 2019.4) pipeline was used to process the sequencing data, as previously described (14). Briefly, raw sequencing reads with exact matches to the barcodes were assigned to the respective samples and identified as valid sequences. The low-quality sequences were filtered through the following criteria (15, 16): sequences that had a length of <150 bp, sequences that had average Phred scores of <20, sequences that contained ambiguous bases, and sequences that contained mononucleotide repeats of >8 bp. Paired-end reads were assembled using FLASH (17). After chimera detection, the remaining high-quality sequences were clustered into amplicon sequence variants (ASVs) with 97% sequence identity using UCLUST. A representative sequence was selected from each ASV using default parameters. ASV taxonomic classification was conducted by BLAST searching the representative sequences set against the Greengenes database (18, 60) using the best hit (19). An ASV table was further generated to record the abundance of each ASV in each sample and the taxonomy of these ASVs. ASVs containing less than 0.001% of the total sequences across all samples were discarded. To minimize the difference in sequencing depth across samples, an averaged, rounded, rarefied ASV table was generated by averaging 100 evenly resampled ASV subsets below 90% of the minimum sequencing depth for further analysis.

Sequence data analyses were mainly performed using QIIME and R (v3.2.0) packages. ASV-level alpha diversity indices, such as the Chao1 richness estimator, abundance-based coverage estimator (ACE) metric, Shannon diversity index, and Simpson index, were calculated using the ASV table in QIIME (20, 61). ASV-level ranked abundance curves were generated to compare the richness and evenness of ASVs among samples. Beta diversity analysis was performed to investigate the structural variation of microbial communities across samples using non-metric multidimensional scaling (NMDS) and the unweighted pair-group method with arithmetic means (UPGMA) hierarchical clustering (21). Differences in the UniFrac distances for pairwise comparisons between groups were determined using Student’s t-test and the Monte Carlo permutation test with 1,000 permutations and visualized through box-and-whiskers plots. Principal coordinate analysis (PCoA) was also conducted based on the genus-level compositional profiles (21). The significance of microbiota structure differentiation between groups was assessed by permutational multivariate analysis of variance (PERMANOVA) (22) an analysis of similarities (ANOSIM) (23, 24) using the R package “vegan.” The taxonomy composition and abundance were visualized using MEGAN (62) and GraPhlAn (25). The R package “VennDiagram” was used to generate Venn diagrams to visualize shared and unique ASVs between samples based on the occurrence of ASVs between samples, regardless of their relative abundance (26). Taxa abundances at the phylum, class, order, family, and genus levels were statistically compared between samples or groups by Metastats (27) and visualized as violin plots. Linear discriminant analysis effect size (LEfSe) was performed to detect differentially abundant taxa between groups using the default parameters (28). Partial least squares discriminant analysis (PLS-DA) was also introduced as a supervised model to reveal the microbiota variation between groups, using the “plsda” function in the R package “mixOmics” (29). Microbial functions were predicted by phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) based on high-quality sequences (30). T-test after inverse sine transformation of ratio data, with a p-value of <0.05 as the criterion for significant difference.

A selective medium (Baird-Parker agar base, Hope Bio-Technology (Qingdao) Co., Ltd.) was used to isolate and culture Staphylococcus from the yolk. The yolks from the two groups were directly coated on the plate medium, placed in a biochemical incubator, and incubated for 48 h at 37°C, respectively. DNA was extracted from the colonies in the culture medium using the Ezup Column Bacteria Genomic DNA Purification Kit (Sangon Biotech Co., Ltd., Shanghai). Sequencing was performed after amplification by PCR, followed by comparative analysis on the NCBI ribosome database (Table 1).

The results showed that the ratio of embryo weight to initial egg weight of group CS was similar to group CC at E7 and E15, but group CS was significantly higher than group CC at E21 (Table 2, p = 0.002, n = 48), and the weight of hatched chicks in group CS was also increased compared to group CC (p = 0.069, n = 48). Fertility (p = 0.549, n = 239), hatchability (p = 0.230, n = 229), and embryonic mortality (p = 0.307, n = 229) did not differ between treatments. To examine the effect of eggshell surface sterilization on the yolk microbiota within fertile eggs, we analyzed the yolk microbial diversity of fertile eggs in the two groups at each sampling time point. The results showed that the alpha diversity of the yolk microbiota in group CS was lower than that of group CC at E07 (Figure 1A) and E15 (Figure 1B), and the difference in all indices between the two groups was significant at E15 (p < 0.050). The PCoA results indicated that the beta diversity of the yolk microbiota in the two groups was similar at the early stages of incubation (Supplementary Figure S1D). However, the clustering of beta diversity of the yolk microbiota was more concentrated in group CS than in group CC at E07 (Figure 1C) and E15 (Figure 1D). By the day the chicken hatched, the diversity of the microbiota in the yolk and gut tended to be the same between the two groups (Supplementary Figures S1E,F).

To analyze the changes in embryo microbiota in the two groups, we compared the relative abundance of the yolk microbiota at different time points, respectively. The results showed that Proteobacteria, Firmicutes, and Bacteroidetes were dominant at the phylum level at all different stages, followed by Actinobacteria and Chloroflexi. The yolk microbiota affected by eggshell surface sterilization is mainly dominated on Proteobacteria and Firmicutes at the phylum level (Figure 2). Notably, the abundance of Proteobacteria in the yolk of the eggshell sterilization group (49.65%) was significantly lower than that of the non-sterilization group (78.53%) at E15 (Supplementary Table S1, p = 0.003, n = 24), while becoming similar at E21 (p = 0.262, n = 24). The abundance of Bacteroidetes in the yolk of eggs with shell sterilization was lower than that of eggs without shell sterilization during embryogenesis (Figure 2), and the same trend was observed in the intestines of hatched chicks (Supplementary Figure S2). At the genus level, Aquabacterium, Acidovorax, Lactobacillus, Novosphingobium, and Azospirillum were the dominant bacteria in the yolk (Figure 2), while Prevotella, Bacteroides, Faecalibacterium, and Clostridium were the dominant bacteria in the intestine (Supplementary Figure S2B). The abundance of Aquabacterium in the yolk of fertile eggs from the surface sterilization group was higher than that of the unsterilized group during the embryogenesis, and the difference reached a significant level at E15 (Supplementary Table S2, p = 0.012, n = 24). The changes in the relative abundance of Acidovorax in the yolk between the two test groups at each developmental stage were similar to those of Aquabacterium. Meanwhile, we found that Lactococcus in the yolk of fertile eggs with surface sterilization was clearly lower than that of the unsterilized eggs at E15 (Supplementary Table S2, p = 0.082, n = 24) and with negligible levels at other time points.

In the previous analysis, we found that the effect of eggshell sterilization before incubation on the yolk microbiota was most pronounced at E15, so we further filtered and analyzed the differential yolk microbes at this time point. We finally filtered out eight major differential microbiota between group CS and group CC, namely, Staphylococcus, Turicibacter, Bdellovibrio, Nitrospira, Microbacterium, Aquabacterium, Pediococcus, and Methylibium (Figures 3, 4). The number of sequences of Staphylococcus (p = 0.009, n = 24), Turicibacter (p = 0.010, n = 24), and Nitrospira (p = 0.012, n = 24) in the yolk of the fertile eggs without surface sterilization was higher than that of sterilized eggs, while the sequence number of Bdellovibrio (p = 0.011, n = 24), Microbacterium (p = 0.031, n = 24), Aquabacterium (p = 0.036, n = 24), Pediococcus (p = 0.125, n = 24), and Methylibium (p = 0.595, n = 24) was lower in the yolk of the fertile eggs without surface sterilization than that in the yolk of sterilized eggs. In addition, LEfSe analysis of the yolk microbiota at different developmental stages showed that the variety of landmark microbiota in the yolk of the unsterilized eggs was higher than that of sterilized eggs at E7, E15, and E21 (Figure S4), especially at E15 (Figure 5). At E15, landmark microbiota at the genus level in the yolk of group CS mainly included Microbacterium, Clostridium, and Asticcacaulis, while that of group CC mainly included Phycicoccus, Arthrobacter, Anaerovorax, Gracilibacter, and Burkholderia (Figure 5).

We performed 16S rDNA full-length sequencing to identify the species of Staphylococcus isolated from the yolk by selective media (Figure S3). The results showed that there were eight species of Staphylococcus saprophyticus and one species of Staphylococcus pseudoxylosus, and the sequence matches were greater than 99.86% (Table 1).