The study included a total of 270 GBS isolates (Supplementary Appendix Figure 1 and Supplementary Appendix Tables 2, 3). In a previous study (Schindler et al., 2020) 240 colonizing isolates and 19 EOD isolates were studied. Of these we randomly chose 126 colonizing, and all EOD isolates for further study. To enrich the sample we also obtained 17 additional EOD isolates from (MHMC) and 108 additional colonizing isolates from Meuhedet Health Maintenance Organization, the 3rd largest health plan in Israel (MHMO). Altogether, we studied 36 GBS isolates obtained from blood cultures of neonates with EOD during 2013–2019, and 234 colonizing GBS isolates collected during routine screening of asymptomatic pregnant women at week 35–37 (Supplementary Appendix Figure 1, Supplementary Appendix Tables 2, 3). The study was approved by the Ethics Committee of MHMC (approval number 0023-18-MHMC).

Group B streptococcus isolates were obtained from clinical cultures and preserved in sterile Brain Heart Infusion broth with 15% glycerol (HY-labs, Israel) at −70°C for long-term storage and at 4°C for short-term maintenance.

Genomic DNA from GBS isolates was extracted using the EZ1 Virus Mini kit (Qiagen) according to the manufacturer instructions. Frozen GBS strains were revived; twice subcultured on trypticase soy with 5% sheep blood agar plates (TSA, HY-labs) at 37°C with 5% CO₂. One to three colonies were suspended in one ml of saline. Following centrifugation at 3,200×g, the bacteria were lysed in a 15 mg/ml lysozyme solution (Lysozyme, Geneaid). A 400 μL of lysed bacteria were transferred for DNA extraction. Extracted DNA was stored at −20°C.

All GBS isolates were serotyped using TaqMan real-time PCR. Amplification was performed using specific primers for every serotype and a fluorescent probe labeled with 6-FAM (Breeding et al., 2016). Positive and negative controls were included in every run.

The presence of genes encoding GBS surface proteins potentially associated with virulence: pili island (PI-1, PI-2a, PI-2b); Rib protein (rib) and hypervirulent GBS adhesin (hvgA) were evaluated among EOD (n = 36), colonizing GBS isolates from MHMC (n = 126) and MHMO (n = 108). The sets of primers used for detection of virulence genes are listed in the Supplementary Appendix Table 3. The reaction mixture, in final volumes of 50 μL, contained: 10 μL of TaqMan fast advanced master mix (Applied biosystems), 0.2 μL of each primer (10 μM), 7.6 μL nuclease-free PCR water, and 2 μL of DNA. Conditions for amplification were as follows: initial denaturation at 94°C for 3 min, followed by 39 cycles of denaturation at 95°C for 30 s, primer annealing at 55°C for 30 s, and extension at 72°C for 30 s. Amplifications were performed in CFX96 thermocycler (Bio-Rad). In each run a negative control consisting of the reaction mixture with nuclease-free water was added. 10 μL of PCR products were visualized by 2% agarose gel electrophoresis with SYBR Safe. The size of each PCR product was estimated by using standard molecular size markers (100-bp ladder) (GeneDireX, Hylabs). The control GBS isolates were used as positive control and was included to each PCR reaction.

A sample of GBS isolates were studied: 24 EOD isolates, and 25 colonizing isolates (Supplementary Appendix Table 4). The genomes of GBS isolates were prepared using Nextera XT kits (Illumina, San Diego, CA, USA) and sequenced using the Illumina MiSeq Reagent Kit v3 (600-cycle). Bioinformatics analysis of GBS whole genome sequencing (WGS) was performed by using the bacterial bioinformatics database and analysis resource PATRIC¹. The reads obtained for each sample were trimmed and the quality of the FASTQ reads was examined using the FASTQ Utilities Service, and finally assembled by SPAdes using the PATRIC website. A high-quality representative genome of S. agalactiae 2603 V/R ATCC BAA611 (serotype V, ST19) was used as reference for genomic comparisons (Tettelin et al., 2002; Ondov et al., 2016). The “Phylogenetic Tree Building” service in PATRIC website provided the Codon Tree method that selects a single copy of the amino acid and nucleotide sequences from a defined number of PATRIC’s global Protein Families (PGFams), picked randomly, to build an alignment, and then generate a tree based on the differences within those selected sequences. We performed the protein sequence-based genome comparisons using bidirectional BLASTP by the PATRIC server. This tool provides information about conserved genomic contexts, and the presence of insertions or deletions. We compared the genomic coding sequences (CDSs) of ST17 strains and ST1 strains to the reference strain (Tettelin et al., 2002). The results of representative ST17 and ST1 strains are displayed with color-coding for protein percent identity relative to the best hit on the reference genome. We searched for virulence factors genes by an exhaustive bioinformatic screening of the database “Virulence Factors Database–VFDB” (Chen et al., 2005), available at the PATRIC website.

Quantitative qRT-PCR analysis of bacterial gene expression was performed as described previously (Sullivan et al., 2017). Primers were designed using Primer3 Plus software and used at a final concentration of 10 μM (Supplementary Appendix Table 5). The GBS isolates used for this study are listed in the Supplementary Appendix Table 4. RNA was extracted from GBS cultures grown at 37°C to an exponential growth phase in BHI medium. RNA was purified using the RNeasy Mini kit (Qiagen) according to manufacturer instructions. Purified RNA was treated with the DNAse kit (HY-labs, Israel) according to manufacturer instructions. cDNA was synthesized using the Hy-RT-PCR kit (HY-labs, Israel), according to manufacturer instructions. cDNA was diluted 1:150 to further reduce bacterial DNA contamination and qPCR was performed using Hy-SYBR power mix (HY-labs, Israel) and CFX96 Real-Time System (Bio-Rad). RNA from three independent biological triplicates were analyzed and final cycle threshold (Ct) for each strain was calculated (the mean value of three experiments). We compared the expression of tested genes between separate groups (EOD vs. colonizing isolates, ST17 vs. ST1 strains, respectively). For this purpose, the final Ct of each group was determined by calculating the mean CT value of all strains belonging to the group. Relative quantification of gene expression was performed using the comparative 2^–ΔΔCT method. The growth rates of colonizing and EOD isolates were tested in BHI medium, and similar growth rates were observed. Results were normalized using rpoB gene as the housekeeping gene. Gene expression results of each gene were normalized to its expression.

Means and medians were computed for continuous variables. Number and percent were used to describe categorical variables. Variables were compared using student T test or chi square as appropriate. P-value was set at < 0.05. For small samples Fisher exact test was used. Mantel–Haenszel common odds ratio estimate was used to calculate 95% confidence intervals.

Multivariate logistic regression (Wald, backward) was performed. The dependent variable was EOD disease. The independent candidate variables were serotypes (3; 6; 1a; 2; 5; 4; other serotypes), and the various virulence factors tested (PI-2a; PI-2b; hvgA; RIB; PI-1).

The sequences are now available to the public at the NCBI (BioProject number BioProject ID PRJNA861829, at the following link: http://www.ncbi.nlm.nih.gov/bioproject/861829).

The dominant GBS serotype isolated from colonized asymptomatic pregnant women in MHMO, was serotype VI [32.4% (35/108 isolates)], followed by serotypes III [22.2% (24/108 isolates)], V [(13.9%, 15/108 isolates)], and IV [10.1% (10/108 isolates)] (Table 1). Serotype II was found only in MHMO isolates. The distribution of GBS serotypes in colonizing isolates in MHMO was comparable to the distribution of GBS serotypes in asymptomatic pregnant women in MHMC (p > 0.05).

We compared the distribution of virulence genes: rib, hvgA, PI-1, PI-2a, and PI-2b among GBS isolates obtained from pregnant women at MHMC (n = 126) and MHMO (n = 108). Overall, the virulence gene distribution among both groups was similar (p > 0.05) (Table 2), except for PI-2a which was less prevalent in MHMC [108/126 (85.7%) vs. 102/108 (94.4%), p = 0.028]. The distribution of virulence genes, among isolates obtained from EOD revealed a high prevalence of surface adhesins compared to isolates obtained from asymptomatic pregnant women (hvgA and rib) (58.3 and 77.8%, respectively; p < 0.01) (Figure 1). Pilus loci PI-2b was also more prevalent among EOD isolates (61.1% vs. 24.4%; p < 0.01), while the pilus loci PI-2a and PI-1 were more frequent among colonizing isolates (89.7 and 93.1% vs. 55.6 and 69.4%; p < 0.01) (Figure 1). The distribution of virulence genes among various serotypes was studied (Table 3). Genes encoding for pilus islands PI-1 and PI-2a were detected in all GBS serotypes, while the gene encoding the pilus loci PI-2b island was mainly detected in serotype III (64.6%). hvgA and rib genes were also detected mainly in serotype III (57.3 and 84.4%, respectively). There was a significant association between hvgA, rib, and PI-2b genes with serotype III (p < 0.01). To determine the frequency of virulence genes among the different STs, we selected 29 EOD and 39 colonizing GBS isolates from MHMC with known sequence types (STs) (32- ST17 isolates, 16 -ST1 isolates, and 20 other ST types) (Supplementary Appendix Table 4). hvgA and rib genes were more prevalent among ST17 isolates (75 and 75.1%, respectively) compared to their distribution in all other groups (p < 0.01). PI-2b island was also more prevalent among ST17 isolates (71.9%, p < 0.01). PI-2a was highly frequent among ST1 isolates (93.7%), while PI-1 island had a similar distribution among all genotypes (Table 4).

We assessed genome sequences of 24 EOD isolates and 25 colonizing GBS isolates from MHMC. General features of the GBS genome sequences, including genome size, number of contigs, and guanine-cytosine (GC) content are summarized in Table 5. We identified a correlation between the ST type and genome size. The genome size of ST17 strains, which was mostly associated with EOD, was significantly smaller 2,001,326 bp compared to the genome size of ST1 strains causing asymptomatic colonization 2,100,416 bp (p < 0.0001).

We constructed a phylogenetic tree of 49 isolates, and a reference strain, to assess the evolutionary relationships among colonizing and EOD GBS isolates. The tree demonstrated five clusters, which corresponded to their MLST STs (Figure 2). The first cluster was mainly composed of ST17, ST19, ST23, and ST27 strains, belonging to EOD and colonizing isolates. The second cluster included mainly colonizing isolates corresponding to various STs (ST1, ST4, ST8, ST12, and ST130). The phylogenetic distribution of the two clusters suggests that they belong to divergent lineages.

To study further the genetic differences between EOD (ST17) and colonizing (ST1) GBS isolates we performed protein sequence-based genome comparisons using bidirectional BLASTP by the PATRIC server. This tool provides information about conserved genomic contexts, and the presence of insertions or deletions. The results of representative nine ST17 and nine ST1 strains are displayed in Figure 3, demonstrate different clusters of conserved genes between the two ST types. There were significant differences in gene sequences between ST17 and ST1 strains. ST1 strains showed a conserved genome relative to the reference strain (>99.9% protein sequence identity), while ST17 strains had less similarity (≤98%). Four mutant gene islands were characterized in both ST17 and ST1 strains (Q1–Q4), compared to the reference strain. In ST17 strains, Q3 and Q4 regions contained both gene deletions and point mutations, while in ST1 strains these regions contained only point mutations. The sequences in these mutant gene islands were mainly coded for phage-related proteins, hypothetical proteins with unknown function, and gene recombination related enzymes. ST17 isolates also carried one mutation-enriched region (Y1), that contained point mutations with sequence identity below 50%. Most of the genes in this region were associated with metabolism and transport (Supplementary Appendix Table 6).

To further characterize the presence and role of virulence genes in GBS isolates, we performed further analyses in 49 sequenced GBS isolates from MHMC by employing the database VFDB using the PATRIC website. All GBS isolates tested shared the following virulence factors, which contribute to adhesion (lmb, hylB), protease activity (scpB), biofilm formation (hasC), and toxins production (cfa, cyl) (Table 6). The main differences between colonizing and EOD isolates were found regarding the Rib surface protein (rib) and fibrinogen binding protein type B (fbsB). These proteins were widespread among EOD isolates and rarely identified among colonizing isolates. Pilus island-2 (PI-2) was frequently identified among colonizing isolates (genes encoded for the backbone pilin protein, PilA and PilC were found in 80% of colonizing isolates). PI-1 had the same distribution across all the isolates regardless of clinical presentation. PI-2 was not detected by this approach.

We performed comparative experiments using quantitative PCR for the expression of the virulence genes for EOD (n = 8) and colonizing (n = 8) isolates. We found that the differences between EOD and colonizing isolates is not only reflected by the presence of virulence genes, but also by their expression (Figure 4). hvgA was barely expressed in colonizing ST17 strains, although the presence of this gene was confirmed by PCR. The expression of the rib gene was two-fold higher in EOD isolates compared to colonizing isolates. Transcription of PI-2a was about three-fold higher in colonizing isolates as compared to EOD isolates. In contrast, the transcription of PI-2b was two-fold higher in EOD compared to colonizing isolates (Figure 5). This suggests difference in the regulation of the pilus loci in different GBS isolates mainly those responsible for invasive (EOD) or non-invasive (colonizing isolates) disease. We also investigated the correlation between the expression levels of PI-2a, PI-2b, and PI-1 genes in ST17 compared to ST1 strains. RT-qPCR analysis was performed with six ST17 and six ST1 GBS strains. Transcription of PI-2a was about two-fold higher in ST1 strains as compared to ST17 strain, while the transcription of PI-2b was two-fold higher in ST17 compared to ST1 strains (Figure 6).

In a multivariate logistic regression analysis PI-I, PI-2A, and serotype 3 were independently associated with EOD (Table 7). The r² was 0.511.