Comparison of the Whole Cell Proteome and Secretome of Epidemic Bordetella pertussis Strains From the 2008–2012 Australian Epidemic Under Sulfate-Modulating Conditions

Sulfate is an important modulator for virulence factor expression in Bordetella pertussis, the causative organism for whooping cough. During infection, sulfate is released when respiratory epithelial cells are damaged which can affect gene expression. The current predominant strains in Australia are found in single nucleotide polymorphism (SNP) cluster I (ptxP3/prn2). It has been reported that ptxP3 strains have higher mRNA expression of virulence genes than ptxP1 strains under intermediate sulfate-modulating conditions (5 mM MgSO4). Our previous proteomic study compared L1423 (cluster I, ptxP3) and L1191 (cluster II, ptxP1) in Thalen–IJssel (THIJS) media without sulfate modulation and identified an upregulation of transport proteins and a downregulation of immunogenic proteins. To determine whether proteomic differences exist between cluster I and cluster II strains in intermediate modulating conditions, this study compared the whole cell proteome and secretome between L1423 and L1191 grown in THIJS media with 5 mM MgSO4 using iTRAQ and high-resolution multiple reaction monitoring (MRM-hr). Two proteins (BP0200 and BP1175) in the whole cell were upregulated in L1423 [fold change (FC) >1.2, false discovery rate (FDR) <0.05]. In the secretome, four proteins from the type III secretion system (T3SS) effectors were downregulated (FC < 0.8, FDR < 0.05) while six proteins, including two adhesins, pertactin (Prn) and tracheal colonization factor A (TcfA), were upregulated which were consistent with our previous proteomic study. The upregulation of Prn and TcfA in SNP cluster I may result in improved adhesion while the downregulation of the T3SS and other immunogenic proteins may reduce immune recognition, which may contribute to the increased fitness of cluster I B. pertussis strains.


INTRODUCTION
Bordetella pertussis is the causative agent for whooping cough, a re-emerging vaccine preventable disease that disproportionately affects infants. The current vaccine used in Australia is the acellular vaccine (ACV) which replaced the previous whole cell vaccine (WCV) in 2000 (Royle and Lambert, 2015). The three-component ACV is the main vaccine used in Australia and contains three virulence factors: pertussis toxin (Ptx), pertactin (Prn), and filamentous hemagglutinin (FHA). In addition to these virulence factors, B. pertussis also produces other virulence factors such as the type III secretion system (T3SS), several toxins [adenylate cyclase toxin (CyaA), tracheal cytotoxin, and dermonecrotic toxin] and several autotransporters involved in adhesion and serum resistance [tracheal colonization factor (TcfA), BipA, serum resistance protein (BrkA), and Vag8] (Melvin et al., 2014).
The expression of these virulence genes is regulated by the BvgAS, a two-component system governing a regulon in response to changes in environmental signals. This regulator responds to changes in environmental stimuli including temperature and the presence of modulators such as sulfate and nicotinic acid (Lacey, 1960;Prugnola et al., 1995;Nakamura et al., 2006). BvgS consists of two Venus flytrap (VFT) periplasmic domains, followed by a Per-Arnt-Sim (PAS) domain and a cytoplasmic histidine kinase (HK) moiety. Two α-helices, linkers 1 and 2, connect the VFT to the PAS domain and HK. BvgS kinase activity is constitutively active under basal conditions when linker 1 adopts a coiledcoil conformation. However, when sulfate or other negative modulators are present, it disrupts the coiled coil conformation of linker 1 which affects the PAS domain interface and causes linker 2 to adopt the coiled coil position resulting in the phosphatase mode (Dupré et al., 2013(Dupré et al., , 2015Lesne et al., 2016Lesne et al., , 2017Lesne et al., , 2018. Under non-modulating (Bvg+) conditions (37 • C and no modulators), BvgS is autophosphorylated. This phosphorylation is then passed through a series of histidine-aspartate transporters and eventually to BvgA. Phosphorylated BvgA will bind to virulence gene promoters and promote transcription (Scarlato et al., 1990). However, when the temperature is 26 • C or modulators (sulfate or nicotinic acid) are present, BvgS and BvgA are not phosphorylated (Bvg−), leading to the repression of virulence genes (Scarlato et al., 1990(Scarlato et al., , 1991. A third phenotype, Bvgi, occurs when there are intermediate modulating conditions which result in lowered levels of phosphorylated BvgA, thereby promoting transcription of a subgroup of virulence genes (Vergara-Irigaray et al., 2005). Different promoters have different affinity for phosphorylated BvgA based on the number of binding sites present. This affects their expression temporally and response to modulators with genes such as fhaB and fim expressed early during Bvg activation while ptx and cyaA are expressed later and are more sensitive to modulation (Scarlato et al., 1991;Kinnear et al., 2001;Veal-Carr and Stibitz, 2005). Bvg+ is necessary and sufficient for infection (Smith et al., 2001), while Bvgi can be important for survival during transmission (Vergara-Irigaray et al., 2005). The role of Bvg− in B. pertussis infection is currently unknown. Bvg− phase-locked mutants were shown to be avirulent (Merkel et al., 1998); however, recent studies suggest that the Bvg− phase may be involved in survival, transmission, and persistence (Hot et al., 2003;Coutte et al., 2016;Moon et al., 2017;van Beek et al., 2018).
Several factors have been linked with the resurgence of pertussis. These factors include improved diagnostic methods, differences in relative vaccine effectiveness, reduced vaccination coverage levels, waning immunity, differences in vaccine immune response, and pathogen adaptation to ACVs (Mooi et al., 2014). It is likely that resurgence is caused by a combination of these factors and the extent of each factor's influence on resurgence is likely to differ in different countries. Despite this, two of the main factors attributed to the increase of pertussis in many countries with high vaccination coverage are the ability of the ACV to protect against disease but not against transmission (Warfel et al., 2014) and pathogen adaptation to the ACV through antigenic mismatches and differences in gene expression (Mooi et al., 2014). It has been shown that ACV immunization elicits predominantly a Th2 response which protects against pertussis disease but does not prevent colonization and transmission of the bacteria to unvaccinated individuals in a baboon model (Ryan et al., 1998;Warfel et al., 2014). Single nucleotide polymorphism (SNP) typing of circulating strains isolated after the introduction of ACVs, found that they primarily belonged to SNP cluster I and contains the alleles ptx promoter 3 (ptxP3) and prn2 (Octavia et al., 2011(Octavia et al., , 2012. SNP cluster I strains have replaced the previously predominant SNP cluster II strains which contained ptxP1/prn3 alleles, and they were responsible for the 2008-2012 Australian epidemic. Notably, ptxP3 strains appeared to produce larger amounts of Ptx (Mooi et al., 2009).
A mixed-infection mice model using representative strains from SNP cluster I (L1423) and II (L1191) showed that cluster I strains were fitter than cluster II in both three-component ACV vaccinated and unvaccinated hosts (Safarchi et al., 2016a). Previous transcriptomic studies have documented differences in the expression of virulence and non-virulence genes between ptxP3 and ptxP1 strains (King et al., 2013). However, our previous comparative global proteomics study between the cluster I and II representative strains showed that there were no differences in the expression of ACV antigens but identified other key protein differences including an increase in TcfA and several transport/binding proteins for phosphate, metals, and amino acids in SNP cluster I and a decrease in immunogenic proteins such as Bsp22, theT3SS tip protein, and lipoprotein BP1569, a toll-like receptor 2 agonist (Luu et al., 2018). Similarly, a study by de Gouw et al. (2014b) also found very little difference in gene expression between Dutch ptxP3 and ptxP1 strains under non-modulating conditions with nine genes significantly downregulated and eight upregulated. However when grown in the presence of intermediate level of sulfate (5 mM MgSO 4 ), ptxP3 strains were found to be less sensitive to sulfate modulation of Bvg and had increased mRNA expression of ACV antigens such as ptx and other virulence factors including vag8 and T3SS. de Gouw et al. (2014b) suggested that during infection, the lysis of host respiratory epithelial cells leads to increased sulfate concentration, possibly through desulfation of host proteins. This results in Bvg modulation and repression of virulence genes (de Gouw et al., 2014b). Decreased sensitivity to modulation may allow B. pertussis to remain virulent for longer periods of time leading to increased transmission and infection (de Gouw et al., 2014b). We hypothesized that differential response to sulfate modulation of Australian SNP cluster I strains compared to SNP cluster II may also contribute to increased fitness observed in vivo (Safarchi et al., 2016a). Therefore, to determine whether additional expression differences exist between the two clusters under modulating conditions, the aim of this study was to use isobaric tag for relative and absolute quantitation (iTRAQ and MRM-hr to compare the global proteome of our representative SNP cluster I and II strain under intermediate sulfate-modulating condition (5 mM MgSO 4 ).

Bacterial Growth Conditions
Bordetella pertussis strains L1423 (cluster I -ptxP3, prn2, ptxA1, and fim3A) and L1191 (cluster II -ptxP1, prn3, ptxA1, and fim3A) from glycerol stocks were grown on Bordet-Gengou (BG, BD scientific) media plates in parallel and incubated at 37 • C for 3-5 days. Pure cultures were obtained by subculturing a single (small, hemolytic) Bvg+ colony onto a second BG plate and incubated again at 37 • C. Pure Bvg+ (small, hemolytic) colonies were then inoculated into THIJS broth (Thalen et al., 1999) with 1% Heptakis [(2,6-O-dimethyl) β-cyclodextrin] and THIJS supplement, and incubated at 37 • C with shaking. After 24 h, the cultures were centrifuged and the cell pellets were resuspended in 100 ml THIJS with 5 mM MgSO 4 and a starting OD 600 of 0.05. The cultures were incubated for 12 h and afterward, the CFU/ml were recorded. Six biological replicates were performed for each strain.

RNA Extraction and RT-qPCR
One milliliter of culture was aliquoted from each sample with 5% v/v phenol in ethanol added to prevent RNA degradation (Bhagwat et al., 2003). RNA was extracted using the ISOLATE II RNA mini kit (Bioline) according to the manufacturer's instructions and further treated with DNase I [New England Biolabs (NEB)] until all traces of genomic DNA were removed. cDNA synthesis was performed using the Tetro cDNA synthesis kit (Bioline). RT-qPCR was performed using the SensiFAST SYBR Hi-ROX kit (Bioline) on the Corbett Rotor-Gene 6000 (Qiagen) to quantify the relative expression of bipA against recA, the housekeeping gene (primers used listed in Supplementary  Table S1) (Stefanelli et al., 2006). The Ct method (Schmittgen and Livak, 2008) was used to calculate the relative expression of bipA between strains grown with and without 5 mM MgSO 4 . Statistical significance was calculated using a two-sided Student's t-test with p < 0.05 defined as significant.
iTRAQ Whole cell and supernatant samples for iTRAQ and MRMhr were prepared as previously described in (Luu et al., 2018; Figure 1). Briefly, centrifugation separated the whole cell and supernatant. The whole cell was then lysed via sonication and the FIGURE 1 | Experimental design of iTRAQ experiment. L1423 and L1191 were grown simultaneously in parallel on BG plates. Pure small, hemolytic Bvg+ colonies were then inoculated into THIJS and then into THIJS + 5 mM MgSO 4 . After 12 h incubation, the whole cell was sonicated and the secretome concentrated. Protein samples were then digested and labeled with iTRAQ reagents. supernatant (100 ml) was concentrated to 200-1000 µl using a 5-kDa Amicon filter (Millipore).
For iTRAQ, two seven-plex iTRAQs were performed for whole cell and supernatant as done previously. Each sevenplex iTRAQ contained three biological replicates from L1423, three biological replicates from L1191, and one pooled reference sample for normalization between iTRAQ experiments. One hundred micrograms of protein from each sample were reduced and alkylated, trypsin digested and labeled with isobaric tags (SCIEX). The labeled samples were then combined, and excess label and salt removed using a cation exchange cartridge followed by an Oasis cartridge clean-up. The combined labeled samples were analyzed using a TripleTOF-5600+ (SCIEX) coupled to a Dionex UltiMate 3000 RSLCnano pump system, Switchos valve unit, and Famos autosampler (Thermo Scientific Dionex). A replicate LC/MS run was performed for each iTRAQ sample.
Protein identification and quantification were performed using ProteinPilot v5.0 (Applied Biosystems) (Shilov et al., 2007) with a custom B. pertussis database. Proteins identifications were accepted with unused score >1.3 (>95% confidence) and >1 peptide. Quantification of proteins was performed using common proteins between both iTRAQs. FCs were calculated in relation to L1423 with L1191 as the reference. Upregulated proteins in L1423 were defined as having FCs >1.2 while downregulated proteins in L1423 had FCs <0.8. Significance was determined using a Student's t-test with FDR correction (Storey and Tibshirani, 2003). p < 0.05 and q < 0.05 was the cutoff for significance. Virulence proteins that showed an upregulated or downregulated (FC > 1.2 or <0.8 and q < 0.1 and q < 0.2) trend were also assessed (Luu et al., 2018).
Cellular location and functional categories were assigned as previously described (Luu et al., 2017). Functional enrichment analysis of proteins was performed using T-profiler adapted to B. pertussis functional categories. T-profiler is a program which was developed for transcriptomic data to determine overall expression differences between predefined genes in functional groups using a T-test (Boorsma et al., 2005). Although originally developed for microarrays, T-profiler has since been applied to proteomics data in Listeria monocytogenes and Escherichia coli O157:H7 (Bowman et al., 2012;Kocharunchitt et al., 2014). T-value and p-value were obtained for each functional category using T-profiler. Briefly, for each predefined protein functional category according to Bart et al. (2010), a two-tailed t-test was performed for each functional category using the mean log2 FC of N proteins assigned to that category against the mean log2 FC of N proteins not assigned to that functional category. To decrease the effect of outliers, the highest and lowest log2 FC value in each gene group were removed prior to calculation as described by Boorsma et al. (2005). Functional categories with T-values >0 were defined as upregulated in cluster I while those with T-values <0 were downregulated (Boorsma et al., 2005). The p-value was corrected for multiple testing using the Benjamini-Hochberg method with p-adjusted < 0.05 assigned as significant (Benjamini and Hochberg, 1995).

MRM-hr
Proteins that were identified to be significantly different in iTRAQ experiments were confirmed using MRM-hr measurements as previously described (Luu et al., 2018). Two or three tryptic peptides from each protein were chosen (Supplementary Table S2) and precursor m/z entered into MRM-hr methods. Each peptide selected was based on the following criteria: 8-20 amino acids in length, no missed cleavage, and no potential ragged ends or contained amino acids susceptible to variable modifications such as cysteine and methionine. "Control proteins" for normalization were also chosen based on the iTRAQ data with additional criteria: FC of 1 ± 0.1 and a coefficient of variation (CV) <10%.
For MRM-hr runs, 10 µg of proteins were digested as described in Luu et al. (2017) and loaded on the LC/TripleTOF 5600+ (SCIEX) and peptides separated over a 30-min gradient. The resulting data were imported into Skyline (v3.5.0.9319) (MacLean et al., 2010); then, intensity measurements from four to six transitions were determined. The inbuilt Skyline R-package, MS-stats (Choi et al., 2014), was then used to determine significantly differentially expressed proteins with linear mixedeffect models (Chang et al., 2012).
To confirm the induction of the Bvgi phenotype with 5 mM MgSO 4 , the relative mRNA expression of bipA was quantified between L1423 and L1191 grown in THIJS with and without 5 mM MgSO 4 . bipA was selected as it is known to be optimally expressed in the Bvgi phenotype (Stockbauer et al., 2001;de Gouw et al., 2014b). The expression of bipA was found to be significantly upregulated when L1191 (FC = 5, p = 0.03) and L1423 (FC = 1.54, p = 0.04) were grown in the presence of 5 mM MgSO 4 , therefore confirming the induction of Bvgi phenotype (Supplementary Figure S1).

Comparison of Whole Cell Proteome Under Intermediate Modulation
From the whole cell, a total of 690 proteins were identified (20% of the total B. pertussis proteome) when B. pertussis was grown under intermediate modulating conditions (5 mM MgSO 4 ). In the first iTRAQ, 603 proteins were identified while in the second iTRAQ, 601 proteins were identified (Supplementary Table S3). There were 513 proteins that were commonly identified, as shown in Figure 2 and Supplementary Table S3, and these proteins were used for quantitative analysis between L1423 and L1191. The functional categories and cellular locations for 690 proteins identified under intermediate modulating condition (5 mM MgSO4) were compared with the 825 proteins previously identified under non-modulating condition (Luu et al., 2018). There was no significant difference in the proportion of functional categories or cellular locations of proteins identified between both conditions (Supplementary Table S4). Furthermore, the proportion of immunogenic (p = 0.08) and Bvgregulated proteins (p = 0.95) also showed no difference between the two conditions with Fisher's exact test.
The expression of proteins in the whole cell between L1423 and L1191 under intermediate modulating conditions were compared using iTRAQ. BP0200, a tripartite tricarboxylate transporter (TCTC) protein, and BP1175, a hypothetical protein with a bacterial oligosaccharide-binding fold (BOF) domain, were found to be significantly upregulated in L1423 (Table 1). Both proteins were also previously identified to be upregulated in L1423 in non-modulating conditions (Luu et al., 2018). No proteins were significantly downregulated. Red circles depict proteins that were significantly upregulated in L1423, and gray circles depict proteins that showed no differences. No proteins were found to be significantly downregulated. # Proteins that were not tested in MRM-hr. $ Proteins that were confirmed to be significantly different using MRM-hr. * Proteins that were unable to be confirmed using MRM-hr. When q value was assessed at 0.1 for virulence proteins that showed a trend in differential expression, Btc22, a T3SS chaperone protein was downregulated in L1423 (FC = 0.75, q = 0.06). At q < 0.2, BcrH1, another T3SS chaperone protein was also downregulated (FC = 0.74, q = 0.13) while BipA showed a trend in upregulation in L1423 compared to L1191 (FC = 1.4, q = 0.13).
For functional enrichment analysis using T-profiler, three functional categories were significantly different. Regulation and ribosome constituents were significantly downregulated in L1423 with a T-value of −5.47 (p-adjust = 0.005) and −4.47 (p-adjust = 0.0005), respectively, while transport/binding proteins were significantly upregulated (T-value = 3.11 and p-adjust = 0.03) (Figure 3A and Supplementary Table S5).

Comparison of Secretome Under Intermediate Modulation
In the secretome, a total ofdownregulation under nonmodulating 173 proteins were detected with 134 proteins identified in the first iTRAQ and 133 proteins in the second iTRAQ. Of these, 96 proteins were commonly identified and used for quantification analysis (Figure 4 and Supplementary  Table S6). Of the 173 proteins, the secretion of 153 proteins could be accounted for bioinformatically and/or were detected in OMVs previously (Roberts et al., 2008;Raeven et al., 2015;Luu et al., 2017). Similar to that seen in the whole cell proteome, 173 proteins identified in 5 mM MgSO 4 were compared to 225 secretome proteins identified in Luu et al. (2018) and there was no difference in the proportion of proteins from different cellular locations and functional categories between cells grown with and without 5 mM MgSO4 (Supplementary Table S7). Furthermore, the proportion of immunogenic proteins were marginally significant [Fisher's exact test (p = 0.052)] but Bvg-regulated proteins [chi square test (p = 0.52)] had no significant differences (Supplementary Table S7).
Under intermediate modulating conditions, 10 proteins were significantly different between L1423 and L1191 ( Table 1). Four proteins (Bsp22, BteA, BopD, and BopN) were significantly downregulated in L1423; all of which belonged to the T3SS. Bsp22 was also previously identified to be significantly downregulated in the secretome of L1423 under non-modulating conditions while BopD and BopN showed a trend of downregulation. Six proteins were significantly upregulated in L1423, two of which were virulence-associated (TcfA and Prn). TcfA was also previously identified to be upregulated in L1423 under non-modulated conditions while Prn was previously borderline significantly upregulated (FC = 1.25, q = 0.053) (Luu et al., 2018). The other four proteins were functionally associated with cell processes (PpiB), energy metabolism (SucD), adaptation (CspA), and cell surface (BP0200). It is interesting to note that BP0200 has been found to be significantly upregulated in both the L1423 whole cell and secretome under non-modulating and intermediate modulating conditions. Expression differences of PtxB and Vag8 between L1423 and L1191 were found to be significant but their FC did not reach the cutoff for downregulation in L1423 at 0.8 (FC = 0.83-0.86 and q = 0.009). When FDR cut-off was analyzed at 0.1 and 0.2, sulphate-binding protein (Sbp) (FC 0.79, q = 0.1) and BipA (FC = 0.74, q = 0.13) were found to show a trend of downregulation in L1423. BipA also showed a trend in downregulation under non-modulating conditions (Luu et al., 2018). There was also an overall downregulation in L1423 of pathogenicity protein category with a T-value of −3.84 (padjust = 0.003) and an upregulation of cell surface proteins with a T-value of 2.77 (p-adjust = 0.02) (Figure 3B and Supplementary  Table S5).

MRM-hr Confirmation of Selected Proteins
For proteins that were identified to be differentially expressed by iTRAQ, MRM-hr was used as a confirmatory technique. Two proteins, BP0200 and BP1175, in the whole cell were increased under intermediate conditions. BP0200 was previously confirmed to be upregulated using MRM-hr in the whole cell and secretome of L1423 under non-modulating conditions while BP1175 was found to be deleted in the genome of L1191 (Luu et al., 2018). Therefore, MRM-hr was only performed for supernatant samples.
Ten proteins were significantly different in the secretome when 5 mM MgSO 4 was added. For four proteins (Prn, PpiB, SucD, and CspA), no suitable peptides for testing were identified. Hence, MRM-hr was only performed on 6 proteins (BteA, BopD, BopN, Bsp22, BP0200, and TcfA) to confirm expression differences. Furthermore, Vag8 and PtxB, which were significant but did not reach the FC cutoff for downregulation in the iTRAQ data were also tested. "Control proteins" for normalization were selected from the iTRAQ data and had the following criteria: FC 1 ± 0.1 and CV < 10%. From these criteria, SucC and Ef-Tu were chosen as control proteins. T3SS proteins including BteA, BopD, BopN, and Bsp22 were all confirmed to be significantly downregulated in L1423 by MRM-hr (FC = 0.09-0.12 and adjusted p-value <0.001) (

DISCUSSION
This study compared the whole cell proteome and secretome of L1423 and L1191, representative strains from SNP cluster I and SNP cluster II, respectively, and identified four proteins downregulated and seven proteins upregulated in L1423 under intermediate sulfate-modulating conditions. Sulfate is an important modulator for virulence factor expression and can be used to induce Bvgi (de Rossi et al., 2001;Friedman et al., 2001;de Gouw et al., 2014a,b,c), and we have confirmed that both of our strains were in the Bvgi state under the 5 mM MgSO 4 modulation condition. The Bvgi intermediate phenotype in B. pertussis has been associated with transmission and biofilm formation (Irie et al., 2004;Vergara-Irigaray et al., 2005;de Gouw et al., 2014c). In this study, we found no significant differences in virulence factor expression in the whole cell proteome but differences in secreted virulence factors were found.
In the whole cell proteome, there were two proteins (BP0200 and BP1175) significantly different between L1423 and L1191 under intermediate modulating conditions. Both proteins were also identified in our previous study under non-modulating conditions to be significantly different in the whole cell, which confirms that regardless of the sulfate conditions, BP0200 and BP1175 are inherently differentially expressed between the two strains (Luu et al., 2018). BP0200 is a member of the tripartite TCTC family with a possible role in nutrient binding (Huvent et al., 2006). Our previous study identified an IS481 element that was present upstream of BP0200 in L1191 but absent in L1423 ( Table 2; Luu et al., 2018). It is possible that IS481 insertion may affect the expression of BP0200. BP1175 contains an oligosaccharide-binding-fold (OB-fold) and is deleted in the genome of L1191 (Luu et al., 2018). Other proteins with this domain has been found in various pathogens, including Salmonella Typhimurium (VisP), Vibrio cholerae (VCA0732), and Pseudomonas aeruginosa (PA0320), to contribute to antimicrobial peptide and/or stress resistance (Fukushima et al., 2012;Moreira et al., 2013;Matson et al., 2017). In S. Typhimurium, VisP inhibits LpxO, an enzyme that modifies lipid A with the addition of a hydroxyl group, and results in increased resistance to antimicrobial peptides such as polymyxin B. Interestingly, B. pertussis also contains a LpxO homolog which suggests that BP1175 may also have a similar role in antimicrobial peptide and stress resistance (MacArthur et al., 2011). Genomic analysis of previously sequenced Australian B. pertussis strains also identified the deletion of BP1175 in 4 other strains belonging to SNP cluster II, V, and unclustered ( Table 2; Bart et al., 2014;Safarchi et al., 2016b). This suggests that the observed difference in BP1175 are due to strain differences rather than cluster difference.
In the secretome, we identified five proteins including two adhesins that were significantly upregulated (TcfA, Prn, CspA, SucD, and PpiB). TcfA was previously found to be upregulated using iTRAQ in the whole cell and secretome of L1423 in unmodulated THIJS broth but could not be previously confirmed using MRM-hr (Luu et al., 2018). The upregulation of TcfA under intermediate modulating conditions was confirmed by MRM-hr ( Table 2). CspA was also previously found under nonmodulating conditions to be upregulated but only in the whole cell. Finally, in our previous study (Luu et al., 2018), Prn was nearly significantly upregulated in the secretome. This study was able to confirm the upregulation of Prn in L1423 under FIGURE 5 | Volcano plot depicting secretome proteins tested using high-resolution multiple reaction monitoring. Red, blue, and gray circles represent upregulated, downregulated, and proteins that were not different, respectively. intermediate modulating conditions. Prn has been shown to be an important adhesin (Van Gent et al., 2011;Safarchi et al., 2015), and it is possible that the upregulation of Prn in L1423 may be associated with differences in Prn type. SNP cluster I strains contain Prn2 while SNP cluster II contain Prn3, and Tohama I, the ACV strain used in Australia contains Prn1 (Octavia et al., 2011). Prn1 contains five GGXXP amino acid repeats in region 1 (R1) while Prn2 contains six repeats in R1 as well as two amino acid differences (Mooi et al., 1998). These differences in the number of repeats were found to affect epitope recognition between antibodies generated against Prn1 and Prn2 (He et al., 2003). However, no difference was observed between Prn3 and Prn1 which contains the same number of amino acid repeats as Prn1 but differs in two amino acids (He et al., 2003). He et al. (2003) proposed that ACV antibodies generated against Prn1 can bind to Prn3 but may be less able to bind to Prn2, thus providing Prn2 strains with a selective advantage to avoid immune recognition. In addition to differences in epitope recognition, Prn1 confers greater colonization ability compared to Prn2 and Prn3, while no difference in colonization ability was observed between Prn2 and Prn3 strains (Van Gent et al., 2011). It is possible that a slight upregulation of Prn combined with Prn2 mismatch with ACV may allow L1423 to compensate for decreased colonization ability of Prn2 in hosts without adversely affecting immune recognition (Safarchi et al., 2016a).
The upregulation of Prn and TcfA, both of which are important adhesins, in our current and previous studies indicates the possibility that SNP cluster I may have increased adhesion and colonization ability compared to SNP cluster II. This hypothesis is supported by the mixed infection mice model results which demonstrated that SNP cluster I were better colonizers than cluster II regardless of the hosts immunization status (Safarchi et al., 2016a).
Several T3SS proteins were found to be significantly downregulated in this study including Bsp22, BopN, BopD, and BteA. Bsp22 forms the T3SS tip complex, while BopD is a translocase for the T3SS effectors BteA and BopN (Fennelly et al., 2008;Nagamatsu et al., 2009;Romero et al., 2013). Bsp22, BopN, and BopD were also previously found to be downregulated in the secretome in non-modulation conditions (Luu et al., 2018). However, previous comparative transcriptomic studies between ptxP3 and ptxP1 strains observed an increase in T3SS transcripts (King et al., 2013;de Gouw et al., 2014b). Moreover, a study comparing non-vaccine and vaccine type isolates in Japan also identified an increase in BteA while a similar study in French isolates identified no BteA expression differences between isolates (Han et al., 2011;Hegerle et al., 2013). Differences in the T3SS result between our study and others may be due to the culture time points used. For example, Han et al. (2011) showed that BteA expression peaks at 24 h, while our proteomic study was conducted at 12 h. It is possible that at an alternative time point, the relative expression of the T3SS between strains may be different. The expression of the T3SS can also be affected by the number of in vitro passages (Gaillard et al., 2011). However, since L1423 and L1191 were cultured in parallel, it is unlikely that reduced T3SS secretion in L1423 is solely due to in vitro passaging.
We examined the T3SS genes for any genetic differences between the two strains used to explain the downregulation of T3SS proteins. No mutations were found within the genes or promoters of the four T3SS genes, bsp22, bopD, bopN, or bteA. In L1423, a non-synonymous SNP was found in bscI at nucleotide position 68 ( Table 2). bscI encodes a T3SS basal body protein homologous to YscI in Yersinia. During the assembly of the T3SS in Yersinia, YscI forms the inner rod protein which allows YscF, the T3SS needle protein, to traverse through the inner membrane (Dewoody et al., 2013). YscI was also found to bind to YscF for assembly of the hollow needle required for the secretion of other T3SS effectors (Cao et al., 2017). The non-synonymous SNP in bscI changed the amino acid from tyrosine to cysteine and was predicted by PROVEAN (v1.1.3) to have a deleterious effect (PROVEAN score = −7.2) on protein function (Choi and Chan, 2015). This suggests that BscI may be impaired in L1423, possibly affecting needle formation for effective secretion of T3SS proteins such as Bsp22, BopD, BopN, and BteA. This hypothesis is supported by studies in Yersinia, where single amino acid mutations to YscI which impaired binding to YscF, was found to reduce the secretion of YopN and LcrV, homologs to BopN and Bsp22, respectively (Cao et al., 2017). Interestingly, genomic analysis has found the same mutation to be present in all SNP cluster I strains as well as other global ptxP3 strains (Bart et al., 2014;Sealey et al., 2014;Safarchi et al., 2016b). Therefore, it is possible that decreased secretion of T3SS effectors may be common across other global ptxP3 strains and may be advantageous in reducing immune recognition to allow current strains to remain hidden from the immune system. Previous in vivo studies have identified Bsp22 and BteA from B. bronchiseptica to be highly immunogenic and protective (Medhekar et al., 2009); however, the homologous proteins in B. pertussis were not identified to be immunogenic or protective in vivo (Hegerle et al., 2013;Romero et al., 2013). Romero et al. (2013) proposed that the inability of T3SS effectors to be protective in vivo may be due to decreased expression or secretion of T3SS effectors in B. pertussis compared to B. bronchiseptica. A western blot by Ahuja et al. (2016) comparing Bsp22 expression in RB50, a B. bronchiseptica strain, with a Tohama I-derived strain and two other Californian clinical B. pertussis strain appears to support decreased expression of Bsp22 in B. pertussis; however, further work is required to confirm this. It is plausible that decreased secretion of T3SS proteins in SNP cluster I may be an ongoing phenomenon to reduce immune recognition in the evolution of B. pertussis from its ancestor, B. bronchiseptica.
In addition to the T3SS proteins, three virulence proteins showed a trend of downregulation (BipA,Vag8,and PtxB) in the supernatant using iTRAQ. Downregulation of PtxB was not confirmed but Vag8 downregulation was confirmed using MRM-hr. A trend of downregulation in BipA was also previously found in the supernatant of L1423 under nonmodulating conditions and was confirmed using MRM-hr (Luu et al., 2018). Both Vag8 and BipA are highly immunogenic proteins with BipA and Vag8 in mice have been shown to be protective in previous immunization studies (de Gouw et al., 2014a,c). Similar to the T3SS, it is possible that the downregulation of BipA and Vag8 may also allow L1423 to lower host immune recognition and contribute to the increased fitness observed in the mouse model (Safarchi et al., 2016a).
Our study did not identify any differences in the expression of sulfatase-containing proteins or sulfate transporters except for a slight but non-significant downregulation of sulfate-binding protein (Sbp) in L1423 (FC = 0.79 and q = 0.1). By contrast, de Gouw et al. (2014b) found that Dutch ptxP3 strains had decreased expression of sulfatase-containing protein genes and increased expression of sulfate transporter genes, such as sbp. It is possible that the differences observed between the two studies may be due to the difference in time points or strains used, or post-transcriptional modifications (Bibova et al., , 2015. Furthermore, the majority of differences seen in this study were observed in the secretome which may not be detectable in transcriptomic experiments if the differences are due to protein secretion and not expression. However, it should be noted that in a recent in vivo transcriptomic study by van Beek et al. (2018), no differences in susceptibility to sulfate modulation response were observed during infection between their ptxP3 and ptxP1 strains, therefore supporting the results in this study.
In this study, we only examined Bvgi condition as 5 mM MgSO 4 induced the Bvgi phenotype in both L1423 and L1191. The Bvg− condition should be investigated in future studies. Recently, RNA sequencing has uncovered hundreds of additional Bvg− genes that are involved in metabolism including sugar and amino acid transport, lipid metabolism, glyoxylate cycle, and phenylacetic acid degradation (Moon et al., 2017). Some of these Bvg− genes were identified to be expressed during in vivo infection, which suggests that the Bvg− phase may play a role in survival, transmission, and persistence (van Beek et al., 2018). Additionally, examining a wider spectrum of modulations such as other MgSO 4 concentrations and additional modulators such as nicotinic acid as well as conditions that mimic the in vivo environment such as low glutamate, iron, and oxygen/carbon dioxide availability may also elucidate more phenotypic differences that contribute to the predominance of SNP cluster I.

CONCLUSION
Whole cell proteome and secretome responses of Australian strains from SNP cluster I (ptxP3) and SNP cluster II (ptxP1) to intermediate levels of sulfate appear to be similar. Although sulfate is an important modulator for B. pertussis virulence gene expression, based on this study, our findings indicate that it is not a major factor for pathogen adaptation and the predominance of SNP cluster I strains in Australia. Of the 11 proteins identified to be significantly different between L1423 and L1191 under intermediate modulating conditions, nine were also differentially expressed in our previous proteomic study under non-modulated conditions (Luu et al., 2018). These include the upregulation of BP0200, CspA, TcfA, and Prn and the downregulation of T3SS and BipA. Notably, the upregulation of Prn and TcfA in SNP cluster I may result in increased colonization ability while the downregulation of the T3SS and other immunogenic proteins may reduce immune recognition, therefore leading to increased fitness of B. pertussis.

AUTHOR CONTRIBUTIONS
LL, SO, and RL designed the study. LL performed the proteomics experiments and drafted the manuscript. LL, SO, LZ, MR, and RL analyzed the results. SO, LZ, MR, VS, and RL provided critical revision of manuscript.

FUNDING
The funding for this project was provided by the National Health and Medical Research Council of Australia (NHMRC).

LL was supported by an Australian Government Research
Training Program Scholarship.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2018.02851/full#supplementary-material FIGURE S1 | RT-qPCR confirmation of Bvgi induction when 5 mM MgSO 4 was added. Relative fold change of bipA in L1191 and L1423 grown with 5 mM MgSO 4 were compared to L1191 and L1423 grown without MgSO 4 , respectively. Error bars displays the 95% confidence intervals, Significant difference with * p < 0.05. FIGURE S2 | Floating bar graph of MRM-hr result for each protein selected. Ef-Tu and SucC were selected as control proteins to normalize against run variation between each sample. The relative abundance of each selected protein in L1191 and L1423 were analyzed with MRM-hr. The thick horizontal line across each bar represents the mean log2 intensity value for the selected protein and was calculated by averaging the mean log2 intensity of each tryptic peptide selected [shown as colored (green, red, and blue) symbols].
TABLE S1 | List of primers used for RT-qPCR.
TABLE S2 | MRM-hr peptides chosen for confirmation of proteins differentially expressed in the secretome between L1423 and L1191.
TABLE S3 | Proteins identified in the whole cell of L1423 and L1191 using iTRAQ. Two sets of seven-plex iTRAQ were performed with three biological replicates in each set. The value shown in columns I-T is the fold change ratio for each biological replicate relative to the pooled reference. The ratios were then averaged for the six biological replicates and divided to obtain the fold change value of L1423 relative to L1191. The accession number, gene and product name, functional category, cellular localization, and whether the proteins are Bvg regulated or previously identified to be immunogenic are also shown.
TABLE S4 | Comparison of the total number of whole cell proteins identified between THIJS and THIJS + 5 mM MgSO 4 in each functional category, cellular localization, immunogenicity, and Bvg regulation. The Fisher's exact or chi square test p-value is also shown.  TABLE S6 | Proteins identified in the L1423 and L1191 secretome using iTRAQ. The value shown in columns K-V is the fold change ratio for each biological replicate relative to the pooled reference. The ratios were then averaged for the six biological replicates and divided to obtain the fold change value of L1423 relative to L1191. The accession number, gene and product name, functional category, cellular localization, and whether the proteins are Bvg regulated, predicted to be secreted, or previously identified to be immunogenic are also shown.
TABLE S7 | Comparison of the total number of secretome proteins identified between THIJS and THIJS + 5 mM MgSO 4 in each functional category, cellular localization, immunogenicity, and Bvg regulation. The Fisher's exact or chi square test p-value is also shown.