ORIGINAL RESEARCH article
Functional Mapping of Phenotypic Plasticity of Staphylococcus aureus Under Vancomycin Pressure
- 1Center for Computational Biology, College of Biological Sciences and Biotechnology, Beijing Forestry University, Beijing, China
- 2College of Biological Sciences and Biotechnology, Beijing Forestry University, Beijing, China
- 3Department of Public Health Sciences and Statistics, Center for Statistical Genetics, The Pennsylvania State University, Hershey, PA, United States
Phenotypic plasticity is the exhibition of various phenotypic traits produced by a single genotype in response to environmental changes, enabling organisms to adapt to environmental changes by maintaining growth and reproduction. Despite its significance in evolutionary studies, we still know little about the genetic control of phenotypic plasticity. In this study, we designed and conducted a genome-wide association study (GWAS) to reveal genetic architecture of how Staphylococcus aureus strains respond to increasing concentrations of vancomycin (0, 2, 4, and 6 μg/mL) in a time course. We implemented functional mapping, a dynamic model for genetic mapping using longitudinal data, to map specific loci that mediate the growth trajectories of abundance of vancomycin-exposed S. aureus strains. 78 significant single nucleotide polymorphisms were identified following analysis of the whole growth and development process, and seven genes might play a pivotal role in governing phenotypic plasticity to the pressure of vancomycin. These seven genes, SAOUHSC_00020 (walR), SAOUHSC_00176, SAOUHSC_00544 (sdrC), SAOUHSC_02998, SAOUHSC_00025, SAOUHSC_00169, and SAOUHSC_02023, were found to help S. aureus regulate antibiotic pressure. Our dynamic gene mapping technique provides a tool for dissecting the phenotypic plasticity mechanisms of S. aureus under vancomycin pressure, emphasizing the feasibility and potential of functional mapping in the study of bacterial phenotypic plasticity.
Phenotypic plasticity is the capacity of an individual genotype to produce different phenotypes (e.g., morphology or behavior) in specific environments (Anderson et al., 2012; Ashander et al., 2016; Kelly, 2019); phenotypic plasticity is a manifestation of biological adaptability and is the main mechanism by which animals and plants respond to some environmental change (Botero et al., 2015; Pozo et al., 2015; Bonamour et al., 2019). A vast number of studies have investigated phenotypic plasticity in plants and animals; however, comparatively few reports have examined the phenotypic plasticity of microorganisms (Savvides et al., 2017). The response mechanisms of microorganisms to environmental change are highly complex. Phenotypic plasticity can mitigate the effects of environmental changes on microbial growth, heredity and development (Moeller and Sanders, 2020). Bacteria show great adaptability to new pressure in changing environments (Sagiv et al., 2015; van Opijnen et al., 2016; Ospina-Serna et al., 2020). Therefore, studying the phenotypic plasticity of microorganisms is key to understanding the response mechanisms of microorganisms to environmental changes.
Genome-wide association studies (GWAS) are an important tool for evaluating phenotypic plasticity. Since the publication of the first successful human GWAS in 2005 (Klein et al., 2005), GWAS have become an increasingly important research method for geneticists. The GWAS method has enabled the characterization of complex diseases and the identification of common genetic variations associated with particular characteristics. However, genetic variation manifests differently in bacteria compared with in humans and other polyploid organisms. Homologous recombination and chromosome segregation occur during human reproduction. On the other hand, bacteria are haploid and asexual with highly structured populations in which horizontal gene transfer and recurrent mutations occur (Read and Massey, 2014). As a result, the development and adaptation of GWAS for use in microorganisms has been relatively slow (Buchanan et al., 2017), although there are some examples of its successful application. GWAS have successfully been applied to identify mutations in the RNA polymerase rpoB gene of Staphylococcus aureus that are associated with drug resistance phenotypes (Alam et al., 2014) and to explore the interactions between Escherichia coli and S. aureus genes (He et al., 2017).
The majority of phenotypic plasticity studies using the GWAS method so far have modeled static phenotypic data measured at a single point in time (Gan et al., 2019). Phenotypes are developmental and changes over time; however, GWAS can only be used to identify phenotypic correlations at a specific time and quantitative trait locus (QTL), without considering dynamic phenotypic development information (Wang et al., 2014; Jiang et al., 2019). Therefore, to assess developmental phenotypic plasticity, dynamic modeling of a series of phenotypes measured over a period of time is required (Li and Wu, 2010). Functional mapping takes into account the biological mechanisms and dynamic processes of phenotypic formation; in this technique, phenotypic traits are continuously measured over a series of time points, and biologically meaningful information such as growth curves is employed to fit these observations. These data are then used to establish a hybrid model under various genetic designs, and the maximum likelihood method and expectation–maximization algorithms are used to estimate parameters. Likelihood ratio statistics have been applied to determine the loci of phenotypic plasticity genes and to quantify their dynamic changes using temporal and spatial scales (Ma et al., 2002; Wu and Lin, 2006). Functional mapping is an accurate and powerful tool for determining the QTL of key growth traits in various species (Li and Sillanpaa, 2015). Bivariate functional mapping is a technique that integrates bivariate data from different conditions into a functional map framework, revealing the expression patterns of phenotypic plasticity genes in different environments (Zhao et al., 2005).
S. aureus is a clinically important bacterial species, which is effectively treated using the antibiotic vancomycin (Shahin et al., 2020; Simon et al., 2020). However, with the continuous use of antibiotics, vancomycin-resistant S. aureus strains are increasingly developing that can adapt to antibiotic environments and show strong phenotypic plasticity (Bakri et al., 2007). In this study, functional mapping and bivariate functional mapping models were applied to evaluate microbial phenotypic plasticity for the first time. We monitored the growth status of S. aureus under four vancomycin concentrations and used these data to establish a model framework with which to study the phenotypic plasticity of microbial development using functional mapping and multivariate analyses. Our model revealed the gene loci and mechanisms related to phenotypic plasticity in S. aureus and for general applications in the study of phenotypic plasticity in microorganisms.
Materials and Methods
Bacterial Strains and Pre-cultivation
First, we measured the vancomycin susceptibilities of 99 vancomycin-sensitive S. aureus strains by subjecting them to vancomycin treatment. Forty-one of the strains were subjected to vancomycin treatment in our previous study (Rong et al., 2019). With the exception of strain S10’, the S. aureus strains showed varying degrees of resistance to vancomycin. Another 58 vancomycin-sensitive strains (S47−S105, with the exception of S86) were clinical isolates from Beijing Chaoyang Hospital and incubated on brain heart infusion (BHI) (Oxoid, Basingstoke, United Kingdom) agar plates supplemented with vancomycin (Sigma-Aldrich, St. Louis, MO, United States) at 50% of the initial minimum inhibitory concentration (MIC) at 37°C. The strains were transferred to fresh medium with the same vancomycin concentration every 24 h for 4 days. Following treatment, we re-calculated the MIC of vancomycin for each strain every 4 days for a total of 60 days, as described previously (Wang et al., 2017). At the end of the treatment, we determined the final MIC value for each strain using the E-test method (bioMérieux); 58 of the strains showed varying degrees of vancomycin resistance. The backgrounds and vancomycin susceptibilities of all strains are listed in Supplementary Table 1. Before the plasticity experiment, single colonies were transferred to BHI medium and shaken overnight at 37°C.
Growth Phenotypic Plasticity Assay
Ninety-nine strains of S. aureus were inoculated into fresh BHI medium supplemented with 0, 2, 4 or 6 μg/mL vancomycin to an initial cellular concentration of 6 × 103 CFU/mL. Growth curve experiments were carried out in 96-well plates (Corning, New York, NY, United States) containing 150 μl volumes in triplicate. Bacterial liquid concentrations were determined by measuring the optical density at 600 nm at 14 time points (1, 2, 4, 6, 8, 10, 12, 16, 20, 24, 30, 36, 42, and 48 h) over a 48 h period. Then, three growth curve equations, Gompertz, Logistic, and Richards, were used to fit the S. aureus developmental phenotype data (Equation 1), and the most suitable equation was selected by least square method. In the following equation, g(T) represents the growth of the strain at time t, A represents the growth degree, R represents the maximum specific growth rate, λ represents the time of lag phase ends.
Whole Genome Sequencing
Genomic DNA was extracted using the TIANamp Bacteria DNA Kit (Tiangen, Beijing, China) according to the manufacturer’s protocol. In total, 99 S. aureus strains were genotyped by genome sequencing on the Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, United States) at Allwegene (Beijing, China). Genome detection, database construction and sequencing, among other factors, may have affected the quality and quantity of the data and therefore affected the results of downstream analyses. To acquire high-quality sequencing data, we performed several quality control steps in the original data, such as removing low-quality reads and joints and calculating the sequencing error rate, Q20 and Q30 statistics, and GC content. S. aureus subspecies NCTC8325 (NC_007795.1) was selected as the control strain for comparison; the comparison results were acquired in BAM format. SAM tools (Li et al., 2009) were used to organize the results, mark repetitive sequences and filter the resulting single nucleotide polymorphisms (SNPs).
Sequencing libraries were constructed using inserted fragment sizes of 300 and 150 bp were sequenced at both ends of the library via paired-end sequencing. The original data size was 697.937200−2658.411600 Mb. After quality control, the data range was 697.16−2581.08 Mb. A total of 110,675 SNPs were obtained from the whole genome data of 99 S. aureus strains following data processing. The sequencing statistics, including the average sequencing data quality and comparisons with the reference sequence, are shown in Supplementary Tables 2, 3. All sequencing data were deposited in NCBI repository1 and the serial numbers are listed in Supplementary Table 1.
Functional Mapping Model
Functional mapping integrates the mathematical relationships between different traits or variables within a genetic mapping framework, where it applies universal biological laws to model the genotypic effects of QTLs. In this study, a functional mapping model was applied to identify key genes affecting the growth of S. aureus at various vancomycin concentrations. Population structure effects were corrected using FastStructure software (Raj et al., 2014), which could be divided into eight groups.
The phenotype of the jth genotype of a QTL was assumed to fall within a multivariate normal distribution:
Here gj was fitted using Equation 1, and the best-fitting equation was selected using the least square method. The fitted mean vector is shown in Equation 3:
The first-order autoregressive (AR) (Yap et al., 2009) was applied to obtain the Σ value. This model comprised the two parameters ρ and σ2:
The likelihood function was then constructed as follows:
Evaluating significant QTL was conducted by hypothesis tests:
H0:(A1,R1,λ1,S1) ≡ (A2,R2,λ2,S2) for j = 1,…,J
H1: Not all equalities above do not hold
Hypothesis test statistics were calculated using the log-likelihood ratio (LR) (Churchill and Doerge, 1994).
where and represent the maximum likelihood estimates under H0 and H1, respectively. The threshold was determined using the permutation test method, in which 1,000 repeated arrangements were made, and the relationship between genotype and bivariate phenotype was randomly reshuffled in each arrangement. A minimum P-value from each arrangement was selected to construct an incremental vector comprising 1,000 elements. Significance was determined at P < 0.05, and the 50th value was taken as the threshold.
Bivariate Functional Mapping Model
A structured antedependence (SAD) model was applied to approximate time correlation covariance matrices of longitudinal traits based on functional mapping. Here, we established a functional mapping method for multiple longitudinal traits. This method can integrate and analyze multiple shapes by applying functional mapping. SAD models are useful tools for designing effective early selection programs for animal and plant breeding, identifying genes that control the progress of human diseases, and raising and solving biological problems at the interface of genetics, development and evolution (Zhao et al., 2005; Wang et al., 2019). The effect of population structure was corrected using FastStructure software (Raj et al., 2014).
The likelihood function of bivariate data of N individuals affected by J genes was formulated by
where w1 and w2 represent the conditional probabilities of genotypes, and fj(zi) represents the corresponding multivariate normal distribution.
Here, gj represents a bivariate mean vector, which was obtained by fitting various processing data from Equation 1.
Then, a first-order SAD (1) model was applied to model Σ.
where Σx and Σy represent covariance matrices of bivariate data, respectively, and Σxy and Σyx represent covariance matrices between two sets of data.
Evaluating significant QTL was conducted by hypothesis tests:
H0: Ω1 ≡ Ω2 for j = 1,…,J
H1: Not all equalities above do not hold
where Ω is the set of all parameters comprising the mean curve. The hypothesis test statistics were calculated via LR, and the threshold was calculated using a permutation test.
Genotype–Phenotype Variation Analysis
Based on the functional positioning and binary curve analyses, the curve parameters fitted by the growth curve equation of different strains were determined. The parameters in this model are as follows: λ Represents the time of the lag phase ends, R represents the maximum specific growth rate of the growth curve, and A represents the growth degree (Jiang et al., 2015). The three growth parameters obtained were biologically significant. The growth curve parameters were compared using t-tests to further explore the significance of the QTLs and how they affect the phenotypic plasticity of S. aureus in relation to growth.
Growth Curve Assay
After 60 days vancomycin treatment, the MICs of most strains were increased. The backgrounds of vancomycin treated strains (S1–S105) and parental strains (Sp1–Sp105) are listed in Supplementary Table 1. All subsequent studies were performed on the vancomycin treated strains, including whole genome sequencing.
The growth curve equation was used to fit the growth data of the 99 S. aureus strains grown in increasing concentrations of vancomycin (0, 2, 4, and 6 μg/mL) and the least square method was used to select the best fitting equation for each individual (Figure 1). The growth curves differed between the strains, which is likely due to their varying MICs and degrees of adaptability to vancomycin treatment. In the absence of antibiotics, all 99 S. aureus strains exhibited typical growth curves (Figure 1A). At a vancomycin concentration of 2 μg/mL, three strains (S22, S25, and S96) were not able to grow. Furthermore, the lag times of 32 strains (S2, S4, S6, S12, S14, S18, S20, S21, S24, S26, S30, S33, S36, S38, S42, S49, S50, S51, S54, S56, S57, S61, S62, S76, S80, S83, S93, S94, S95, S97, S104, and S105) were significantly delayed, and the time to reach the logarithmic growth phase was longer than that without antibiotic supplementation (Figure 1B). At a vancomycin concentration of 4 μg/mL, 18 strains (S1, S2, S4, S6, S12, S14, S17, S20, S22, S24, S25, S33, S36, S42, S49, S54, S95, and S96) did not grow, and the lag times and stationary phases of 35 strains were delayed (Figure 1C). Thirty-two strains (S1, S2, S4, S5, S6, S7, S12, S14, S17, S19, S20, S21, S22, S24, S25, S26, S29, S30, S31, S33, S36, S42, S49, S54, S62, S77, S88, S91, S94, S95, S96, and S103) did not grow at 6 μg/mL group due to the high vancomycin concentration. For each strain, differences between four growth curves under different vancomycin concentrations could directly reflect growth phenotypic plasticity (Figure 1). After equation fitting, three key parameters (A, R, and λ) were collected (Supplementary Table 4). The curve equations could associate with gene data and locate the SNPs that affect the dynamic growth process by functional mapping.
Figure 1. Growth curve fitting of 99 S. aureus strains. The growth curves of all strains at 0 μg/mL (A), 2 μg/mL (B), 4 μg/mL (C), and 6 μg/mL (D) vancomycin treatment were fitted by Gompertz, Logistic or Richards equation. Each line represents the growth of a strain, red represents Gompertz equation, blue represents Richards equation, and green represents Logistic equation.
To locate candidate genes related to S. aureus resistance development and phenotypic plasticity at different vancomycin concentrations, the dynamic growth data of 99 S. aureus strains were analyzed by functional mapping. Thirty-eight significant SNPs were identified under the four vancomycin concentrations (Figure 2), and their gene locations were functionally annotated, including 16 missense mutations and 16 synonymous mutations (Table 1).
Figure 2. Manhattan plots of functional mapping. Manhattan plot of the significance test based on log-likelihood ratio (LR) for associations between the phenotypic growth differences in different groups and SNPs distributed on the S. aureus genome. SNPs positions in genome were divided by two colors. Each dot in the plot represented an SNP, and a reference line was used on the y-axis to reflect genome-wide significance. (A–D) Represent Manhattan plots of four treatment groups.
SNP 25261 is located within the SAOUHSC_00020 gene, which is related to the transcriptional regulatory protein WalR. SAOUHSC_00169 (containing SNP 183485) encodes the ABC transporter permease, which is involved in the transport of various molecules such as amino acids, lipids, lipopolysaccharides, inorganic ions, peptides and sugars; this transporter may play a crucial role in bacterial drug resistance. SNP 193712 is located within SAOUHSC_00176, a chemoreceptor that recognizes constituents of transport systems and initiates signal transduction pathways. SAOUHSC_00544 (harboring SNP 550323) encodes fibrinogen-binding protein SdrC. SNP 2728755 is located in SAOUHSC_02967, which is involved in the formation of biofilms, arginine synthesis and pH regulation. SAOUHSC_00199 (containing SNP 220967) encodes acyl-CoA, which also plays an important role in the formation of S. aureus biofilms. SNP 1961488 is located in SAOUHSC_02078, which encodes phi PV83 orf 10-like protein. SNP 2771775 is located in SAOUHSC_02998, which is involved in the synthesis of the outer capsule of the cell wall. SAOUHSC_00069 (containing SNP 73870) encodes immunoglobulin G-binding protein A, which plays an important role in the inhibition of the host innate and adaptive immune responses. SAOUHSC_02026 (containing SNP 1930634), SAOUHSC_02029 (containing SNP 19304358), SAOUHSC_02030 (containing SNP 1935188), and SAOUHSC_02060 (containing SNP 1953238) encode phi ETA orf 58-like protein, phi ETA orf 56-like protein, phi ETA orf 55-like protein and phi PVL orf 51-like protein, respectively. The SAOUHSC_01121 gene (containing SNP 1076587) is related to alpha-hemolysin, which binds to the membrane of eukaryotic cells, resulting in the release of low-molecular weight molecules, eventually leading to osmotic lysis. Alpha-hemolysin also inhibits host neutrophil chemotaxis to the bacterial lesion. SAOUHSC_02182 (containing SNP 2045672) encodes the tail tape measure protein. The SNPs at positions 229384, 2072412, 448470, and 942441 were located within non-coding regions.
Bivariate Functional Mapping
Functional mapping focuses on the genes that affect the dynamic growth of strains in different concentrations of antibiotics using univariate phenotype data. Growth curves in the presence or absence of antibiotics formed the typical bivariate phenotype, while bivariate functional mapping can simultaneously analyze correlated phenotypes in plasticity study. The differences in the growth curves of S. aureus cultured with and without vancomycin directly reflected the phenotypic plasticity of S. aureus. Therefore, we performed bivariate functional mapping to analyze binary data collected following S. aureus growth at vancomycin concentrations of 0 and 2 μg/mL, 0 and 4 μg/mL, and 0 and 6 μg/mL (Figure 3). A total of 40 significant SNP sites were found, including 6 missense mutations and 26 synonymous mutations (Table 2). The SNP at position 2062927 is located within the SAOUHSC_02214 gene, which encodes a conserved hypothetical phage protein. SAOUHSC_00025 (containing SNP at site 31826) encodes an uncharacterized GRAM_POS_ANCHORING domain-containing protein. SNP 183439 is located in SAOUHSC_00169, which encodes the same permease transporter protein discussed above. The SNP at position 1928590 is located in SAOUHSC_02023, an important player in cell separation, daughter cell formation and autolysis. The SAOUHSC_02033 gene (containing SNP 1938367) is related to the phage tape measure protein. The SNP at position 1942424 is located in SAOUHSC_02036, which encodes a phage structural protein. The SAOUHSC_02871 gene (containing SNP 2644158) encodes an acetyltransferase.
Figure 3. Manhattan plots of bivariate functional mapping. Manhattan plot of the significance test based on the associations between phenotypic plasticity of microbial growth in the environment with or without antibiotics and SNPs. Each dot in the plot represented an SNP, and a reference line was used on the y-axis to reflect genome-wide significance. (A–C) Bivariate analysis of Manhattan plots with 0 and 2 μg/mL, 0 and 4 μg/mL, or 0 and 6 μg/mL vancomycin treatment.
Genotype–Phenotype Variation Analysis of Growth Parameters vs. Significant SNPs
For the genotype–phenotype variation analysis, we selected seven key genes affecting the phenotypic plasticity of S. aureus under various vancomycin concentrations based on their P-values and occurrence frequencies from functional mapping and bivariate functional mapping. The growth curves of two genotypes corresponding to the seven selected genes were fitted (Figure 4), and the parameters of the curve were assessed by t-tests to explore the specific effects of the genes on phenotypic plasticity (Table 3). We found that there were significant differences in the growth curves of the representative genotypes.
Figure 4. Growth variations of significant SNPs. For each SNP, strains were divided into two groups according to their genotypes. The thick red and blue lines represent the average growth curves of the two genotypes, while the thin lines in the background represent the individual growth curves.
The SAOUHSC_00020 (SNP25261) gene was divided into two groups according to base A and G, and the growth of strains with genotype A was significantly higher than that of individuals with genotype G; the difference in growth rates between the two genotypes increased with the vancomycin concentration. A similar pattern was observed for SAOUHSC_00544 and SAOUHSC_00176. For the SAOUHSC_02998, SAOUHSC_00025, SAOUHSC_00169, and SAOUHSC_02023 genes, significant differences in growth between the two representative genotypes were present only at high vancomycin concentrations. We therefore categorized the first gene group as global regulator plasticity genes and the second group as stress response plasticity genes.
The SAOUHSC_00020 gene (containing SNP 25261) significant in three concentrations (0, 4, and 6 μg/mL) of function mapping and is therefore likely to represent a key gene affecting S. aureus growth. The degree of growth and maximum specific growth rate were significantly different between the two SAOUHSC_00020 genotypes examined, while there was no significant difference in the lag phase length. SAOUHSC_00176 (containing SNP 193712), which encodes a protein involved in extracellular solute binding, was significant in two sets of functional mapping (0 and 4 μg/mL); all growth parameters were significantly affected by SAOUHSC_00176. SAOUHSC_00544 (containing SNP 550323) encodes the fibrinogen-binding protein SdrC, which plays an important role in adhesion and pathogenesis. SAOUHSC_00176 also affected the entire growth process in all directions. The SAOUHSC_02998 gene (containing SNP 2771777) encodes the capsular polysaccharide biosynthesis protein Cap5C and can trigger an increase in growth rate. The SAOUHSC_00025 gene (containing SNP 31826) encodes GRAM_POS_ANCHORING domain-containing protein, SAOUHSC_00169 (containing SNP 183485) encodes an ABC transporter permease, and SAOUHSC_02023 encodes a bifunctional autolysin; these three genes are key regulators of phenotypic plasticity, as indicated by significant differences in all growth parameters at a vancomycin concentration of 6 μg/mL.
Phenotypic plasticity enables microorganisms to grow and develop in different environments, playing an important role in the competition and evolution of microorganisms. However, studies investigating phenotypic plasticity in microorganisms are limited (Zheng et al., 2020). In our previous work, we used a bivariate GWAS method to locate genes that affected the phenotypic plasticity related to the growth of S. aureus at a single time point (Rong et al., 2019). However, this was not sufficient to evaluate the entire complex bacterial growth process. In this study, we analyzed the phenotypic plasticity of 99 S. aureus strains grown under different antibiotic concentrations using functional mapping and bivariate functional mapping. We identified seven genes that were significantly related to growth-related phenotypic plasticity under the pressure of vancomycin in S. aureus.
With the increasing use of antibiotics, more drug-resistant S. aureus strains are emerging, posing a serious threat to human public health (Laupland, 2013). Vancomycin resistant S. aureus (VRSA) can adapt to antibiotic environments (Gardete and Tomasz, 2014) via cell wall thickening, reducing intracellular toxicity, reducing peptidoglycan cross-linking and modifying autolysis rates. Previous studies indicated that spontaneous mutations could significantly affect the evolution of vancomycin intermediate resistant S. aureus (VISA) (Chen et al., 2014; Hu et al., 2016). Phenotypic plasticity plays an important role in bacterial adaptive evolution, which highlights the importance of investigating the mechanisms driving phenotypic plasticity in S. aureus during vancomycin exposure.
So far, the majority of GWAS studies have explored the relationships of genes with bacterial growth and development by analyzing a single phenotype. Applying bivariate GWAS improved the scope of this method and reduced the rate of false positives (Bedo et al., 2014); however, both GWAS and bivariate GWAS are performed using a single time point. As bacterial development is a dynamic process, bacterial phenotypic plasticity studies should take into consideration the whole growth process. GWAS cannot be used to fully and accurately determine the relationship between genes and phenotypes throughout the development process (Wang et al., 2014). Our functional mapping technique based on dynamic time series data showed great potential for future applications in plasticity research.
In this study, 99 S. aureus strains were cultured under varying vancomycin concentrations, and 78 significant loci were identified by functional mapping and bivariate functional mapping. Seven genes that play significant roles in S. aureus phenotypic plasticity were identified: SAOUHSC_00020(walR), SAOUHSC_00176, SAOUHSC_00544(sdrC), SAOUHSC_02998, SAOUHSC_00025, SAOUHSC_00169, and SAOUHSC_02023. WalR plays a role in the development of bacterial drug resistance (Baseri et al., 2021; Zhu et al., 2021) and regulates genes involved in cell wall metabolism, virulence regulation, biofilm production, oxidative stress resistance and antibiotic resistance via the direct or indirect regulation of autolysins (Kan et al., 2007; Delaune et al., 2012; Poupel et al., 2016). Single nucleotide substitutions within the walR gene lead to vancomycin and daptomycin co-resistance and cause the typical cell wall thickening observed in antibiotic-resistant clinical isolates (Howden et al., 2011). SAOUHSC_00176 is involved in ATP transport and cell membrane synthesis, playing a crucial role in S. aureus colonization and growth (Oiki et al., 2019). SdrC is a cell surface-associated calcium-binding protein that regulates adhesion and pathogenesis and can promote bacterial adhesion by mediating interactions with components of the extracellular matrix (Batool et al., 2021; Wang et al., 2021). The SAOUHSC_02998 gene (containing SNP 2771777) encodes a protein involved in the biosynthesis of capsular polysaccharides and plays a key role in resisting adverse environments and enhancing cellular adhesion (Liu et al., 2018). SAOUHSC_00169 (SNP 183485) encodes ABC transporter permease, which is involved in the transportation of ATP and in drug resistance. An important mechanism driving bacterial antibiotic resistance is the reduction of intracellular antibiotic concentrations by multi-drug resistance proteins (Baltes et al., 2020; Gao et al., 2020). SAOUHSC_02023 (SNP 1928590) encodes bifunctional autolysin, which regulates cell division, daughter cell separation, and autolysis (Albrecht et al., 2012; Kluj et al., 2018). SAOUHSC_00025 (SNP 31826) encodes an uncharacterized GRAM_POS_ANCHORING domain-containing protein that was identified to play an important role in phenotypic plasticity regulation by our bivariate functional mapping model. The correlation between MIC and growth inhibition of each strain under vancomycin pressure was also analyzed, but not statistically significant (P = 0.12) (data not shown). These results indicated dynamic gene mapping techniques are more accurate and powerful.
Our genotype–phenotype growth variation analysis enabled categorization of the seven significant SNPs into two groups: global regulator genes and stress response plasticity genes. For the global regulator genes (SAOUHSC_00020, SAOUHSC_00544, and SAOUHSC_00176), one genotype exhibited superior growth traits under all vancomycin concentrations (Figure 4). For the stress response plasticity genes (SAOUHSC_02998, SAOUHSC_00025, SAOUHSC_00169, and SAOUHSC_02023), the difference in growth between the two genotypes was significant only under high vancomycin concentrations, at which one phenotype exhibited a significant growth advantage and improved plasticity (Figure 4). At the same time, the association between genotypes and MICs of these seven genes was analyzed. The result showed the genotypes of three genes (walR, SAOUHSC_00176, and sdrC) were associated with MICs, while the other four had no significant correlation (Supplementary Figure 1).
In this study, we analyzed the phenotypic plasticity of S. aureus under vancomycin pressure by functional mapping, which led to the identification of seven significant plasticity genes. Functional mapping and bivariate functional mapping provide novel strategies for the study of bacterial phenotypic plasticity, enabling analysis of the relationships between multiple phenotypes and genotypes throughout the dynamic development process. Our method facilitated the location of significant genes, which would not have been possible using a standard, static GWAS methodology.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The name of the repository can be found below and the accession numbers can be found in the Supplementary Material: NCBI repository, https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA722566.
YJ and RW conceived and designed the experiments. DY and XZ performed the experiments. LJ, MY, and XH analyzed the data. DY and YJ wrote the manuscript. All authors reviewed the manuscript.
This work was supported by the Fundamental Research Funds for the Central Universities (2021ZY64), the Beijing Municipal Funds for Talent Training (2017000020124G276), and the Natural Science Foundation of China (31971398).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.696730/full#supplementary-material
Alam, M. T., Petit, R. A., Crispell, E. K., Thornton, T. A., Conneely, K. N., Jiang, Y., et al. (2014). Dissecting vancomycin-intermediate resistance in Staphylococcus aureus using Genome-Wide Association. Genome Biol. Evol. 6, 1174–1185. doi: 10.1093/gbe/evu092
Albrecht, T., Raue, S., Rosenstein, R., Nieselt, K., and Goetz, F. (2012). Phylogeny of the Staphylococcal major autolysin and its use in genus and species typing. J. Bacteriol. 194, 2630–2636. doi: 10.1128/jb.06609-11
Anderson, J. T., Inouye, D. W., Mckinney, A. M., Colautti, R. I., and Mitchell-Olds, T. (2012). Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change. Proc. R. Soc. Lond. Ser. B-Biol. Sci. 279, 3843–3852. doi: 10.1098/rspb.2012.1051
Ashander, J., Chevin, L.-M., and Baskett, M. L. (2016). Predicting evolutionary rescue via evolving plasticity in stochastic environments. Proc. R. Soc. Lond. Ser. B-Biol. Sci. 283:20161690. doi: 10.1098/rspb.2016.1690
Bakri, F. G., Abu Al-Hommos, N., Shehabi, A., Naffa, R. G., Cui, L., and Hiramatsu, K. (2007). Persistent bacteraemia due to methicillin-resistant Staphylococcus aureus with reduced susceptibility to vancomycin in a patient with erythrodermic psoriasis. Scand. J. Infect. Dis. 39, 457–460. doi: 10.1080/00365540601105780
Baltes, F., Pfeifer, V., Silbermann, K., Caspers, J., Von Rekowski, K. W., Schlesinger, M., et al. (2020). beta(1)-Integrin binding to collagen type 1 transmits breast cancer cells into chemoresistance by activating ABC efflux transporters. Biochim. Biophys. Acta Mol. Cell Res. 1867:118663. doi: 10.1016/j.bbamcr.2020.118663
Baseri, N., Najar-Peerayeh, S., and Bakhshi, B. (2021). The effect of subinhibitory concentration of chlorhexidine on the evolution of vancomycin-intermediate Staphylococcus aureus and the induction of mutations in walKR and vraTSR systems. Infect. Genet. Evol. 87:104628. doi: 10.1016/j.meegid.2020.104628
Batool, N., Shamim, A., Chaurasia, A. K., and Kim, K. K. (2021). Genome-Wide Analysis of Staphylococcus aureus sequence type 72 isolates provides insights into resistance against antimicrobial agents and virulence potential. Front. Microbiol. 11:613800. doi: 10.3389/fmicb.2020.613800
Bonamour, S., Chevin, L.-M., Charmantier, A., and Teplitsky, C. (2019). Phenotypic plasticity in response to climate change: the importance of cue variation. Philos. Trans. R. Soc. B-Biol. Sci. 374:20180178. doi: 10.1098/rstb.2018.0178
Botero, C. A., Weissing, F. J., Wright, J., and Rubenstein, D. R. (2015). Evolutionary tipping points in the capacity to adapt to environmental change. Proc. Natl. Acad. Sci. U.S.A. 112, 184–189. doi: 10.1073/pnas.1408589111
Buchanan, C. J., Webb, A. L., Mutschall, S. K., Kruczkiewicz, P., Barker, D. O. R., Hetman, B. M., et al. (2017). A Genome-Wide Association Study to identify diagnostic markers for human pathogenic campylobacter jejuni strains. Front. Microbiol. 8:1224. doi: 10.3389/fmicb.2017.01224
Chen, C. J., Lin, M. H., Shu, J. C., and Lu, J. J. (2014). Reduced susceptibility to vancomycin in isogenic Staphylococcus aureus strains of sequence type 59: tracking evolution and identifying mutations by whole-genome sequencing. J. Antimicrob. Chemother. 69, 349–354. doi: 10.1093/jac/dkt395
Delaune, A., Dubrac, S., Blanchet, C., Poupel, O., Maeder, U., Hiron, A., et al. (2012). The WalKR system controls major Staphylococcal virulence genes and Is involved in triggering the host inflammatory response. Infect. Immun. 80, 3438–3453. doi: 10.1128/iai.00195-12
Gao, Q., Li, X.-X., Xu, Y.-M., Zhang, J.-Z., Rong, S.-D., Qin, Y.-Q., et al. (2020). IRE1 alpha-targeting downregulates ABC transporters and overcomes drug resistance of colon cancer cells. Cancer Lett. 476, 67–74. doi: 10.1016/j.canlet.2020.02.007
He, X., Jin, Y., Ye, M., Chen, N., Zhu, J., Wang, J., et al. (2017). Bacterial genetic architecture of ecological interactions in co-culture by GWAS-taking Escherichia coli and Staphylococcus aureus as an example. Front. Microbiol. 8:2332. doi: 10.3389/fmicb.2017.02332
Howden, B. P., Mcevoy, C. R. E., Allen, D. L., Chua, K., Gao, W., Harrison, P. F., et al. (2011). Evolution of multidrug resistance during Staphylococcus aureus infection involves mutation of the essential two component regulator WalKR. PLoS Pathog. 7:e1002359. doi: 10.1371/journal.ppat.1002359
Hu, Q. W., Peng, H. G., and Rao, X. C. (2016). Molecular events for promotion of vancomycin resistance in vancomycin intermediate staphylococcus aureus. Front. Microbiol. 7:1601. doi: 10.3389/fmicb.2016.01601
Jiang, L., Clavijo, J. A., Sun, L., Zhu, X., Bhakta, M. S., Gezan, S. A., et al. (2015). Plastic expression of heterochrony quantitative trait loci (hQTLs) for leaf growth in the common bean (Phaseolus vulgaris). New Phytol. 207, 872–882. doi: 10.1111/nph.13386
Jiang, L., Sung, L., Ye, M., Wang, J., Wang, Y., Bogard, M., et al. (2019). Functional mapping of N deficiency-induced response in wheat yield-component traits by implementing high-throughput phenotyping. Plant J. 97, 1105–1119. doi: 10.1111/tpj.14186
Kan, J., Suzuki, M. T., Wang, K., Evans, S. E., and Chen, F. (2007). High temporal but low spatial heterogeneity of bacterioplankton in the Chesapeake bay. Appl. Environ. Microbiol. 73, 6776–6789. doi: 10.1128/aem.00541-07
Klein, R. J., Zeiss, C., Chew, E. Y., Tsai, J. Y., Sackler, R. S., Haynes, C., et al. (2005). Complement factor H polymorphism in age-related macular degeneration. Science 308, 385–389. doi: 10.1126/science.1109557
Kluj, R. M., Ebner, P., Adamek, M., Ziemert, N., Mayer, C., and Borisova, M. (2018). Recovery of the peptidoglycan turnover product released by the autolysin Atl in Staphylococcus aureus Involves the Phosphotransferase System Transporter MurP and the Novel 6-phospho-N-acetylmuramidase MupG. Front. Microbiol. 9:2725. doi: 10.3389/fmicb.2018.02725
Liu, J., Yang, L., Hou, Y., Soteyome, T., Zeng, B., Su, J., et al. (2018). Transcriptomics study on Staphylococcus aureus biofilm under low concentration of ampicillin. Front. Microbiol. 9:2413. doi: 10.3389/fmicb.2018.02413
Ma, C.-X., Casella, G., and Wu, R. (2002). Functional mapping of quantitative trait loci underlying the character process: a theoretical framework. Genetics 161, 1751–1762. doi: 10.1093/genetics/161.4.1751
Oiki, S., Sato, M., Mikami, B., Murata, K., and Hashimoto, W. (2019). Substrate recognition by bacterial solute-binding protein is responsible for import of extracellular hyaluronan and chondroitin sulfate from the animal host. Biosci. Biotechnol. Biochem. 83, 1946–1954. doi: 10.1080/09168451.2019.1630250
Ospina-Serna, J., Huber, P., Odriozola, M., Fermani, P., and Unrein, F. (2020). Picocyanobacteria aggregation as a response to predation pressure: direct contact is not necessary. FEMS Microbiol. Ecol. 96:fiaa153. doi: 10.1093/femsec/fiaa153
Poupel, O., Moyat, M., Groizeleau, J., Antunes, L. C. S., Gribaldo, S., Msadek, T., et al. (2016). Transcriptional analysis and subcellular protein localization reveal specific features of the essential WalKR system in Staphylococcus aureus. PLoS One 11:e0151449. doi: 10.1371/journal.pone.0151449
Pozo, M. J., Lopez-Raez, J. A., Azcon-Aguilar, C., and Garcia-Garrido, J. M. (2015). Phytohormones as integrators of environmental signals in the regulation of mycorrhizal symbioses. New Phytol. 205, 1431–1436. doi: 10.1111/nph.13252
Read, T. D., and Massey, R. C. (2014). Characterizing the genetic basis of bacterial phenotypes using genome-wide association studies: a new direction for bacteriology. Genome Med. 6:109. doi: 10.1186/s13073-014-0109-z
Rong, M., Zheng, X., Ye, M., Bai, J., Xie, X., Jin, Y., et al. (2019). Phenotypic plasticity of Staphylococcus aureus in liquid medium containing vancomycin. Front. Microbiol. 10:809. doi: 10.3389/fmicb.2019.00809
Sagiv, J. Y., Michaeli, J., Assi, S., Mishalian, I., Kisos, H., Levy, L., et al. (2015). Phenotypic diversity and plasticity in circulating neutrophil subpopulations in cancer. Cell Rep. 10, 562–573. doi: 10.1016/j.celrep.2014.12.039
Savvides, A., Van Ieperen, W., Dieleman, J. A., and Marcelis, L. F. M. (2017). Phenotypic plasticity to altered apical bud temperature in Cucumis sativus: more leaves-smaller leaves and vice versa. Plant Cell Environ. 40, 69–79. doi: 10.1111/pce.12835
Shahin, I. G., Abutaleb, N. S., Alhashimi, M., Kassab, A. E., Mohamed, K. O., Taher, A. T., et al. (2020). Evaluation of N-phenyl-2-aminothiazoles for treatment of multi-drug resistant and intracellular Staphylococcus aureus infections. Eur. J. Med. Chem. 202:112497. doi: 10.1016/j.ejmech.2020.112497
Simon, A., Moreira, M. L. A., Costa, I. F. D. J. B., De Sousa, V. P., Rodrigues, C. R., Da Rocha, E., et al. (2020). Vancomycin-loaded nanoparticles against vancomycin intermediate and methicillin resistantStaphylococcus aureus strains. Nanotechnology 31:375101. doi: 10.1088/1361-6528/ab97d7
van Opijnen, T., Dedrick, S., and Bento, J. (2016). Strain dependent genetic networks for antibiotic-sensitivity in a bacterial pathogen with a large pan-genome. PLoS Pathog. 12:e1005869. doi: 10.1371/journal.ppat.1005869
Wang, J., Zhang, M., Wang, M., Zang, J., Zhang, X., and Hang, T. (2021). Structural insights into the intermolecular interaction of the adhesin SdrC in the pathogenicity of Staphylococcus aureus. Acta Crystallogr. F-Struct. Biol. Cryst. Commun. 77, 47–53. doi: 10.1107/s2053230x21000741
Wang, Q., Gan, J., Wei, K., Berceli, S. A., Gragnoli, C., and Wu, R. (2019). A unified mapping framework of multifaceted pharmacodynamic responses to hypertension interventions. Drug Discov. Today 24, 883–889. doi: 10.1016/j.drudis.2019.01.009
Wang, Y., Li, X., Jiang, L., Han, W., Xie, X., Jin, Y., et al. (2017). Novel mutation sites in the development of vancomycin-intermediate resistance in Staphylococcus aureus. Front. Microbiol. 7:2163. doi: 10.3389/fmicb.2016.02163
Yap, J. S., Fan, J., and Wu, R. (2009). Nonparametric modeling of longitudinal covariance structure in functional mapping of quantitative trait loci. Biometrics 65, 1068–1077. doi: 10.1111/j.1541-0420.2009.01222.x
Zhao, W., Hou, W., Littell, R. C., and Wu, R. (2005). Structured antedependence models for functional mapping of multiple longitudinal traits. Stat. Appl. Genet. Mol. Biol. 4:33. doi: 10.2202/1544-6115.1136
Zheng, X., Bai, J., Ye, M., Liu, Y., Jin, Y., and He, X. (2020). Bivariate genome-wide association study of the growth plasticity of Staphylococcus aureus in coculture with Escherichia coli. Appl. Microbiol. Biotechnol. 104, 5437–5447. doi: 10.1007/s00253-020-10636-6
Keywords: phenotypic plasticity, Staphylococcus aureus, functional mapping, vancomycin, growth curve
Citation: Yang D, Zheng X, Jiang L, Ye M, He X, Jin Y and Wu R (2021) Functional Mapping of Phenotypic Plasticity of Staphylococcus aureus Under Vancomycin Pressure. Front. Microbiol. 12:696730. doi: 10.3389/fmicb.2021.696730
Received: 17 April 2021; Accepted: 23 August 2021;
Published: 09 September 2021.
Edited by:David W. Ussery, University of Arkansas for Medical Sciences, United States
Reviewed by:Hanne Ingmer, University of Copenhagen, Denmark
Zulfiqar Mirani, Pakistan Council of Scientific and Industrial Research, Pakistan
Copyright © 2021 Yang, Zheng, Jiang, Ye, He, Jin and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yi Jin, firstname.lastname@example.org