Species-Wide Phylogenomics of the Staphylococcus aureus Agr Operon Revealed Convergent Evolution of Frameshift Mutations

Abstract

Staphylococcus aureus is a prominent nosocomial pathogen that causes several life-threatening diseases, such as pneumonia and bacteremia. S. aureus modulates the expression of its arsenal of virulence factors through sensing and integrating responses to environmental signals. The agr (accessory gene regulator) quorum sensing (QS) system is a major regulator of virulence phenotypes in S. aureus. There are four agr specificity groups each with a different autoinducer peptide sequence encoded by the agrD gene. Although agr is critical for the expression of many toxins, paradoxically, S. aureus strains often have nonfunctional agr activity due to loss-of-function mutations in the four-gene agr operon. To understand patterns in agr variability across S. aureus, we undertook a species-wide genomic investigation. We developed a software tool (AgrVATE; https://github.com/VishnuRaghuram94/AgrVATE) for typing and detecting frameshift mutations in the agr operon. In an analysis of over 40,000 S. aureus genomes, we showed a close association between agr type and S. aureus clonal complex. We also found a strong linkage between agrBDC alleles (encoding the peptidase, autoinducing peptide itself, and peptide sensor, respectively) but not agrA (encoding the response regulator). More than 5% of the genomes were found to have frameshift mutations in the agr operon. While 52% of these frameshifts occurred only once in the entire species, we observed cases where the recurring mutations evolved convergently across different clonal lineages with no evidence of long-term phylogenetic transmission, suggesting that strains with agr frameshifts were evolutionarily short-lived. Overall, genomic analysis of agr operon suggests evolution through multiple processes with functional consequences that are not fully understood. IMPORTANCE Staphylococcus aureus is a globally pervasive pathogen that produces a plethora of toxic molecules that can harm host immune cells. Production of these toxins is mainly controlled by an active agr quorum-sensing system, which senses and responds to bacterial cell density. However, there are many reports of S. aureus strains with genetic changes leading to impaired agr activity that are often found during chronic bloodstream infections and may be associated with increased disease severity. We developed an open-source software called AgrVATE to type agr systems and identify mutations. We used AgrVATE for a species-wide genomic survey of S. aureus, finding that more than 5% of strains in the public database had nonfunctional agr systems. We also provided new insights into the evolution of these genetic mutations in the agr system. Overall, this study contributes to our understanding of a common but relatively understudied means of virulence regulation in S. aureus.

Keywords: Staphylococcus aureus; agr; bioinformatics; convergent evolution; genomics; hemolysis; mutation rate; phylogenetics; quorum sensing.

FIG 1

Distribution of agr groups across 40,890 S. aureus genomes from the Staphopia database. AgrVATE was used to assign the agr groups and genomes with unknown agr groups were filtered out. (A) A schematic depiction of the agr operon showing two divergent promoters (P2 and P3) driving agrBDCA and small RNA RNAiii. (B) Frequency of each agr group in the Staphopia database. (C) Frequency of agr groups across the major clonal complexes (CC) of S. aureus. (D) Relative proportions of agr groups from S. aureus isolated from different body sites in percentage.

FIG 2

AgrA evolves independently of agr group with only two major amino acid sequence configurations across S. aureus. (A) Scoring agr group exclusivity in clusters of unique AgrABC amino acid sequences. Extracted amino acid sequences of each agr gene were clustered with 100% identity to obtain all possible AA sequence configurations. Each cluster was then scored based on the number of agr groups the cluster sequence was found in (1 = one agr group; 4 = four agr groups) represented by a circle. Only clusters with more than 50 sequences are shown. The color of each circle represents the number of sequences within the cluster. The red and blue arrows indicate the major (AgrA^K136) and minor clusters (AgrA^R136) of AgrA AA sequences, respectively. (B) Amino acid sequence alignment of the two major alleles of AgrA. (C) Maximum likelihood phylogeny (GTR+FO model, 1000 ultrafast bootstrap replicates with average bootstrap support of 97.8%) of 334 S. aureus strains with each tip representing a unique ST. Tip colors represent the AgrA alleles and the corresponding heatmaps show the agr group and clonal complex of each tip. Scale bar indicates the number of substitutions per site. All tips representing the AgrA^R136 allele are confined to the clade highlighted in blue. (D) Linkage disequilibrium (LD) block plot of the agr operon and 1000 bp flanking regions. Each point on the block indicates R² values of LD calculated by plink for a given pair of SNPs. The y-axis indicates the distance between SNP pairs.

FIG 3

Presence of putative non-functional variants of the agr operon. (A) Frequency of frameshift mutations in coding regions of unique agr operon sequences across the Staphopia database. Arrows indicate agr genes and bars indicate the number of frameshifts at the corresponding position (bin width = 40). Bar colors represent each agr group. (B) Frequency (right) and effect (left) of commonly occurring frameshift mutations across unique agr operon sequences. Bars are colored based on agr group and arrows are colored based on agr gene, black outlines represent canonical protein length and red outlines represent truncated protein lengths. Labels (center) indicate the amino acid change due to the frameshift mutation. (C) Normalized percentage of samples with noncanonical two-component regulator (TCS) gene lengths. Histidine kinase (HK) and response regulator (RR) genes of TCS were extracted from the Staphopia database and commonly occurring gene lengths (>5000 genomes) were excluded. The remaining strains were considered to have noncanonical gene lengths.

FIG 4

Identical agr mutations evolve independently of phylogeny across different clonal complexes. (A) Bars show the frequency of mutations at a given site in descending order. Frequency of each mutation in a dereplicated set of CC8, CC30, CC22, and CC5 samples (left). Frequency of each mutation in a set of randomly selected CC8/CC30/CC22/CC5 agr operons with simulated indels (right). (B) The minimum number of changes on the tree versus the number of occurrences of frameshift mutations. Each circle represents a position on the agr operon that has acquired a mutation in at least 2 samples in a dereplicated set of CC8, CC30, CC22, and CC5 sequences (left). The x-axis represents the number of times the position has acquired a frameshift mutation. The consistency index and the minimum number of changes on the tree were measured at these sites for each CC from the respective phylogenetic tree (GTR+FO model, 1000 ultrafast bootstrap replicates with average bootstrap support of at least 71%) Blue line follows y = x distribution. The outlier CC8 point (black arrow) corresponds to a previously characterized agrA mutation that was not a true agr null (59). The consistency index and the minimum number of changes on the tree were measured for phylogenies for the respective CCs where the tree tips were randomly shuffled (right). One hundred shuffled trees were generated per CC.