compGWAS: a new GWAS tool allows revelation of the genetic architecture and risk stratification for the versatile pathogen Streptococcus pyogenes

doi:10.21203/rs.3.rs-5049272/v1

Download PDF

Research Article

compGWAS: a new GWAS tool allows revelation of the genetic architecture and risk stratification for the versatile pathogen Streptococcus pyogenes

https://doi.org/10.21203/rs.3.rs-5049272/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Gene inactivations caused by loss-of-function mutations and regulatory changes caused by insertions/deletions (InDels) are common genetic factors closely related to phenotypic diversity or pathogenic versatility of many bacterial species. However, these genetic factors were usually ignored by the computational approach of genome-wide association studies (GWAS). It prevents the full understanding of the contributions of genetic variants to phenotypic diversity or the roles in shaping genetic architecture of bacterial species of diverse phenotypes. Group A Streptococcus pyogenes (GAS) is one of the most versatile pathogens causing a variety of primary diseases, as well as disease progressions, complications, and sequelae and is a promising species to do investigations in this regard.

Methods

By using GAS as a paradigm, we developed a new GWAS tool, compGWAS, to comprehensively identify phenotype-associated genetic variants that include not only SNPs, but also InDels and gene inactivations. The genetic architecture of GAS phenotypes was revealed by considering all these types of variants. A GWAS polygenic score (GPS) model was developed through integration of all types of associated variants for phenotype stratification.

Results

By leveraging this newly developed tool, we constructed a relationship network between 1,361 variants linked with 783 genes and eight GAS phenotypes. The network shows a high level of polygenicity of the GAS phenotypes (ranging from 6 to 148 genes) and pleiotropicity of the causal genes (as many as eight phenotypes). Further investigation revealed a unique genetic architecture of GAS phenotypes as a combination of many low-effect common variants and a small proportion of high-effect low-frequency variants with gene inactivations being predominant. By adding gene inactivations and InDels, the proportion of explained phenotypic variance increased by 7%-16%, resulting in a total explained variance as high as 50%. The high explained variance allowed us to construct a GPS model with high discriminatory capabilities in GAS phenotype stratification with the AUC > 80% in the validation dataset.

Conclusions

Our work provides a novel tool and analysis framework for investigating phenotypic effects and genetic characteristics of InDels and gene inactivations previously ignored. Our study has implications for understanding genetic architecture of versatile pathogens like GAS.

Group A Streptococcus pyogenes

Genome-wide association study (GWAS)

Genetic architecture

Phenotype prediction

Group A Streptococcus pyogenes (GAS) is a leading global pathogen exclusively causing human diseases [1]. Its disease manifestations are highly diverse, ranging from superficial skin or pharyngeal infections to more invasive and severe tissue/organ failures [2]. The severe diseases caused by GAS include, but are not limited to, rheumatic fever, necrotizing fasciitis, streptococcal toxic shock syndrome, and sepsis. A peculiar property of GAS infection is that many cases of the infections are symptomatic at early onset, can can become serious very quickly and cause lethal outcomes if not diagnosed or treated in a proper manner and time-frame [3, 4]. A deep understanding of the genetic architecture of infection phenotypes could guide early diagnosis and clinical prevention of GAS disease progression. A delineation of the genetic differences between GAS phenotypes is a critical component of their genetic architecture.

Previous studies have identified a series of genetic factors related with GAS pathogenesis, such as the cell-surface adhesion M protein, the covRS regulatory system, the cysteine protease (speB), the streptolysin O (slo), and the IL8 peptidase (scpC) [2, 5, 6], the correlation is largely strain-specific or emm type-specific. However, there is still a large gap between GAS genomics and pathogenic phenotypes. The systematic connections between genetic factors and GAS phenotypes remain to be defined. Many questions regarding genetic architecture of the complex GAS phenotypes are still unrevealed, such as which combination of the genotypes may affect specific phenotypes, the effect sizes of the genotypes, and the contributions of different genomic regions to these effects.

Genome-wide association studies (GWAS) are promising approaches to delineate these questions [7]. Previous studies employed GWAS to GAS genomes of specific emm types [8] or specific phenotypes [9, 10] to identify genetic factors associated with GAS pathogenesis. The major limitation of these studies, and many other microbial GWAS studies, is that they essentially used single variant-based GWAS (or SNP-based GWAS), thus only allowing tests for associations of single nucleotide variants with the phenotypes of interest [11, 12]. This limitation hinders the comprehensive identification of genotype-phenotype associations in GAS and other pathogenic bacterial species, and further undermines the deep understanding of the genetic architecture of the pathogens.

Single variant-based GWAS analysis is insufficient to detect other complex forms of variants underlying the phenotypic diversity and pathogenic complications for many microbial genomes, including GAS, which are enriched with genetic variability. Firstly, this approach lacks the power to detect the associations of genetic loci that have diverse allelic forms but lead to convergent functional outcomes. One notable example is the virulence regulator gene, covS, encoded by GAS. Its inactivation has been experimentally confirmed to be highly associated with hypervirulence of GAS strains and this inactivation can be caused by a variety of Insertion/Deletion (InDel) mutations in a strain-dependent manner, such as a single bp deletion of T in a serotype M1 isolate, MGAS5005, and M53 isolate, AP53 [13, 14], a five bp deletion of GAAAA in a serotype M23 isolate M23ND [15], an eleven bp insertion in a serotype M81 isolate, UH328 [16]. However, the single variant-based GWAS analyses when applied to GAS genomes were unable to detect this association between covS inactivation and hypervirulence of GAS [8–10]. Another example is the InDels in segments of low sequence complexity, which are prone to alignment ambiguities. The alignment ambiguities can introduce false isomorphs from identical alleles and consequently lose the power to detect the association for the true allele. Secondly, the single variant-based GWAS fails to consider the synergistic effect of multiple mutations on specific genetic units, such as the gene-coding sequences (CDS) or non-coding sequences (nonCDS), which might be more directly linked to functional alterations than single mutations. One illustrative scenario is the mutations in intergenic regulatory regions, where multiple genetic changes at distinct loci in a specific region may together produce regulatory effects. Novel GWAS methods are needed to overcome these issues by considering the complex variant forms to inform more associations on the whole gene level or functional unit level.

Several k-mer-based GWAS methods were proposed to address the issues by testing for associations between k-mers of variable lengths and phenotypes [17–20]. This approach automatically accounts for multiple mutations on the k-mers. However, the convergent outcomes of the k-mer mutations on genetic units cannot be directly evaluated and the biological interpretation of the k-mer testing results is still a challenging task.

In the present study, we developed a new GWAS analysis tool containing two new methods, i.e., IndelWinScan and CDS2GWAS, which allow for association tests for genetic units of CDS and nonCDS that harbor complex variant forms. We offer the new GWAS tool in an open-source package, compGWAS, that can be broadly applied to other microbial GWAS analyses. By integrating the new GWAS methods with the traditional single variant-based GWAS, we performed comprehensive identifications of genotype-phenotype associations for a set of GAS phenotypes. For the first time we characterized a polygenic genetic architecture of GAS phenotypes in a functionally interpreted manner. More importantly, by leveraging the polygenic nature of GAS phenotypes, we established composite polygenic score models for GAS phenotype stratification in a machine learning framework yielding high predicting performance.

Genomic data preparation and preprocessing

The genomic data of 2598 GAS strains were collected from the NCBI Genbank database (ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank). The demographic and clinical information of the strains (i.e. isolation country, and isolation year, isolation tissue, and clinical phenotype) was obtained via extensive search against public databases and existing literatures. The main database sources include NCBI (https://www.ncbi.nlm.nih.gov/data-hub/genome/?taxon=1314), ATCC (www.atcc.org), and NCTC (https://www.culturecollections.org.uk/collections/nctc.aspx). A large proportion of the strains derived their clinical information from other work [10].

Two subsets of GAS were selected from the whole set of 2598 strains as the training and validation dataset, respectively. Specifically, 790 strains were selected as the training cohort with balanced emm serotype distribution and full-range coverage of isolation time points, and 710 strains were selected as the validation cohort with similar emm serotype distribution and time coverage to that of the training cohort.

Emm serotyping

The emm serotyping of the strains were performed by BLAST against the M protein database from CDC (ftp://ftp.cdc.gov/pub/infectious_diseases/biotech) and a set of representative Mga regulon proteins we cataloged for emm pattern A-E. The emm type of each strain was determined based on the sequence similarity with the reference M proteins and the genomic organization of the Mga regulon proteins. Furthermore, the phylogenetic tree was constructed based on M protein sequences to confirm the serotyping results by assuming that the strains on close branches tend to have the same emm type in comparison to the strains on distant branches. The spurious serotyping was re-examined and corrected based on the phylogenetic relationship of the strains.

Data cleaning, core genome generation, and mutation detection

The genomes were filtered based on the following criteria: (i) the assembly genome is shorter than 90% of the median length (1,646,681 bp); (ii) the genome assembly is too fragmented (the number of contigs longer than 500 bp > 50); (iii) large fragments other than phage elements were missing; (iv) the number of missing or interrupted genes > 15; (v) the emm serotype was unable to be unambiguously determined. Finally, 136 genomes were excluded and 2466 were retained for further analysis. Among them, 165 strains with complete genome and high assembly quality were used for generating the core genome.

The core genome was created by aligning the GAS genomes against one of the complete genomes, AP53, which was sequenced and annotated with high quality by us previously [21]. The core genomic sites were defined as those covered by more than 98% of the 165 selected complete genomes. The highly variable regions were excluded, including phage elements, transposons, and repetitive regions.

The mutations were detected using varScan (version 2.3.8) [22] and the functional annotation of the mutations was performed using the annotation module in the package “compGWAS” we developed. The phylogenetic tree was constructed based on the mutations in core genome and visualized using iTOL [23].

Genotype transformation for the new GWAS methods CDS2GWAS and IndelWinScan

In order to perform association tests for the coding status of CDS and mutational status of non-coding region using our newly developed methods CDS2GWAS and IndelWinScan, we first transformed the coding status of CDS and mutational status of non-coding regions to categorical genotypes.

The coding status of each CDS was transformed to binary phenotypes of intact:inactivated. Initially, initial gene locations on GAS genomes were derived from orthologous sequence alignments with the reference genome. Then the sequences in the aligned regions along with 165bp upstream and downstream of the aligned regions were extracted from each genome. Finally, the coding status of the extracted sequences was determined based on the following rules: (i) the CDS encodes an intact complete gene if the CDS has a length of multiple of 3 and contains a start codon (ATG/TTG/CTG/TTG) at the alignment start site and stop codon (TAA/TGA/TAG) in the alignment end site without internal stop site; (ii) the CDS encodes an inactivated gene if the CDS has a length of multiple of 3 and contains a start codon at the beginning and stop codon in the end, as well as an internal codon; (iii) the CDS encodes an intact complete gene if the CDS lacks a start codon at the beginning or a stop codon at the end, but contains a start codon within ± 12bp from the starting site or a stop codon within the ± 24bp from the end site, and the length between the start codon and end codon is multiple of 3.

The mutational status of non-coding regions was transformed to binary phenotypes of intact:mutated induced from InDels. The InDels were called and the InDel minor allele was identified at each alignment loci. A sliding window strategy was used to detect the occurrence of any InDel minor allele in a window of fixed length with a certain step size. Within each window, if the InDel minor alleles occurred in one or more loci for a strain, the mutational status of this fragment in this strain was considered to be mutated, otherwise is intact. The window size was taken as 200bp and step size as 20bp. If the window size is longer than the length of the intergenic region, then the window only covered the intergenic region.

GWAS analysis

Three modalities of GWAS were performed, i.e., the conventional SNP-based GWAS for SNPs, the newly developed CDS-based GWAS for the transformed genotypes of coding status of CDS, and InDel-based GWAS for the transformed genotypes of mutational status of nonCDS. The genotypic variants in the form of SNPs, coding status of CDS, and mutational status of nonCDS were filtered with minor allele count < 2 and allele number > 2. The highly divergent genes were also filtered (number of mutation sites/gene length > 15%). The association test was performed by fitting a binary logistic regression model:

$$\:{Y}_{logit}=ln\left(\frac{P(Y=\text{0,1})}{1-P(Y=\text{0,1})}\right)={Z}_{\beta\:}\left(X\right)={\beta\:}_{0}+{\beta\:}_{1}{X}_{1}+{\beta\:}_{2}{X}_{2}+{\beta\:}_{3}{X}_{3}$$

where Y represents the binary phenotypes, X₁ represents the biallelic genotypes. β₁ is the regression coefficient of the genotypes X₁. β₀ is the intercept. The logistic model allows adjusting for the effects of confounding factors. Here, X₂ and X₃ are the confounding factors population structure and GAS collection time, respectively. β₂ and β₃ are the regression coefficients of the two confounding factors, respectively. The test significance of the regression coefficients was evaluated using Wald test.

Determination of the association p-value thresholds

The significance p-value thresholds for SNP-based GWAS in each comparison group were determined separately using the maximal curvature method. Basically, the linear curve of p-values in descending order was fitted using the polynomial model and the maximal curvature point (or the elbow point) was used as the cutoff for the association p-value [24]. Considering the low number of variants from CDS-based and InDel-based GWAS, the p-value thresholds for the two categories were determined manually based on Manhattan plots.

Variant screening for genomic partitioning and functional analysis

LD-based pruning and thresholding strategy was used to remove the redundancy induced from correlated variants for each GWAS analysis. LD analysis was performed using Haploview [25] and the LD blocks across the genome were detected for SNPs with pair-wise correlation r² ≥ 0.8 within a distance of 10 kb (Figure S6). Then, the SNPs were screened in a hierarchical way. For the coding regions in each block, only the non-synonymous SNPs with the strongest association p-value were retained for positively associated and negatively associated SNPs, respectively. For the non-coding regions in each block, the SNPs with the strongest association p-value were retained for positively associated and negatively associated SNPs, respectively. For the SNPs not belonging to any LD block, only non-synonymous SNPs were retained. The genes harboring the selected non-synonymous SNPs were used for further analysis.

For InDel mutational fragments from InDel-based GWAS, only the fragments with the strongest association p-value in each non-coding intergenic region were used. The genes potentially regulated by the InDel mutational fragments were defined as the downstream genes with the 5’-ends within 1000 bp from the center of the fragments (5’-nonCDS) or the upstream genes with the 3’-ends within 500 bp from the center of the fragments (3’-nonCDS) were used for functional analysis. The distance of 1000 bp and 500 bp here was chosen based on the common lengths of promoter regions and termination regions for prokaryotes.

The truncated CDS from CDS-based GWAS with association p-value lower than the general thresholds described in the previous section were directly used for functional analysis.

Construction of network diagrams for the association analysis resultsThe network relationship between the genotypes and phenotypes were constructed on the gene level for three gene sets: (i) the genes harboring non-synonymous SNPs but without being truncated; (ii) the genes potentially regulated by the non-coding SNPs and InDel mutational fragments located in the 5’-end (5’-nonCDS) or 3’-end (3’-nonCDS) of the genes; (iii) the genes truncated by SNPs or InDels within the genes. The positive and negative association relationships were separately considered and indicated in edge color. The weights of the network edges indicate the association significance. For the genes in set (i), the weights are proportional to –log10P of the most significant non-synonymous SNPs. For the genes in set (ii), the weights are proportional to the –log10P of the most significant InDel mutational fragment or non-coding SNPs. For the genes in set (iii), the weights are proportional to the –log log10P of the inactivated gene. The network relationship between the functional categories and phenotypes were constructed by integrating the effects of multiple genes within each pathway on the phenotype. The combined association p-value p_comb of multiple genes within each pathway was evaluated using Fisher’s product test with the assumption that the association tests for individual genes in the same functional category are independent [26, 27]:$\:{p}_{comb}={\prod\:}_{i=1}^{n}{p}_{i},\:\:$and $\:F=-2*{\sum\:}_{i=1}^{n}\text{l}\text{n}\left({p}_{i}\right)$follows the $\:{\chi\:}_{2n}^{2}$ distribution. The weights of the network edges are proportional to -log₁₀(p_comb). The functional category node sizes are proportional to the number of genes in each category. The positive and negative association relationships were separately considered and indicated in edge color. All network diagrams were generated using Cytoscape [28].

Estimation of variance explained by the genetic variants

Least Absolute Shrinkage and Selection Operator (LASSO) regression in the R package glmnet (version 2.0.16) was employed to select the variants without co-linearity. The LASSO regression model between GAS phenotypes and multiple genetic variants was fitted using L₁ regularization:

$$\:{Y}_{logit}=ln\left(\frac{P(Y=\text{0,1})}{1-P(Y=\text{0,1})}\right)={\theta\:}_{0}+{\sum\:}_{j=1}^{n}{\theta\:}_{j}{g}_{j}$$

where Y represents the binary phenotype, g_j is the j-th variant, θ_j represents the regression coefficient for the j^th variant. θ₀ is the intercept. The variants with regression coefficient smaller than 10^− 5 were removed and the remaining variants were used for calculation of the explained variance and AUC.

The proportion of variance explained was estimated using the metric “coefficients of determinant” (R²) on the liability scale [29]. We first fitted the penalized logistic regression model between GAS phenotypes and variants with the R package logistf [30]:

which is able to reduce the bias of maximum likelihood estimates induced from complete separation (a predictor variable separates outcome variables completely without the need for estimation of the model parameters). Then the R² was calculated on the logit liability scale [31]:

$$\:R2\:=\:\text{v}\text{a}\text{r}\left({Y}_{logit}\right)\:/\:\left(\text{v}\text{a}\text{r}\right({Y}_{logit})+{\pi\:}^{2}/3)$$

where var(Y_logit) is the variance of Y on the logit liability scale.

Construction of GPS prediction models and performance evaluation

The composite GPS was developed as a weighted sum of the carried minor allele across all variants for individual samples. The weight is defined by the regression coefficient of the corresponding variant:

$$\:\text{G}\text{P}\text{S}={\sum\:}_{i}^{n}allele\left(i\right)\times\:coeff\left(i\right),$$

where n is the number of genetic variants used for the prediction model, allele(i) represents the minor allele carrier status of the i^th variant for the sample (allele = 1 for carrier, and allele=-1 for non-carrier), and coeff (i) is the regression coefficient for the i^th variant. The final GPS was fitted to a regular logistic regression model with the phenotype as the outcome in the training cohort of GAS:

$$\:ln\left(\frac{P(Y=\text{0,1})}{1-P(Y=\text{0,1})}\right)=PRS$$

The GPS model was then used to predict the phenotypes of GAS strains in the validation cohort. The prediction performance of the model was evaluated using the ROC analysis to calculate the AUC values with the R package “pROC”.

Dataset overview

We established a large-scale GAS genomics dataset of 2,598 GAS strains from public databases and published literature (Additional File 1). The dataset encompasses more than 25 pathogenic phenotypes and 152 different emm serotypes from 29 isolation countries. The clinical information and geographical distribution of the top 14 phenotypes of the GAS strains was presented in Fig. 1. Three notable features were observed: (i) no single phenotype was dominated by strains of a single emm serotype (Fig. 1a), (ii) the majority of the strains were collected in the last two decades when high-throughput sequencing technologies became available (Fig. 1b), and (iii) the strains are broadly distributed across global geographic locations (Fig. 1c). Those features make the current GAS genomics data suitable for systematic analysis of statistical association between phenotypes and genotypes.

We identified the core genome of GAS using a selected set of 165 strains with complete genomes and high sequence quality (Additional File 2), yielding 1,566,429 core genomic sites distributed in 1584 core genes. The majority of the core genes have a high level of mapping coverage among the GAS strains and 138 genes have a low level of coverage (the number of strains whose gene coverage less than 50% > 7) (Figure S1). We excluded 104 of the 138 genes from association tests due to the low mapping coverage (Additional File 3).

Mutation detection

We selected a training cohort of 790 strains with balanced emm serotype distributions and full-range coverages of isolation time points for mutation detection (Additional File 4) and focused on the top eight phenotypes containing a sufficient number of strains, i.e., invasive (INV), pharyngitis (PHY), superficial soft tissue infection (SSTI), superficial soft tissue infection (SSTI), deep soft tissue infection (DSTI), rheumatic fever (RF), sepsis (SP), scarlet fever (SF), and pneumonia (PM), for model fitting and association tests (Fig. 1d-1e). The phylogenetic structure based on the GAS core genome shows that the strains do not cluster by phenotypes, emm serotypes, or isolation year, and therefore the population stratification or epidemiological factors will not significantly influence the statistic association test (Figure S2). Furthermore, our method assesses the relative allele frequencies only between paired phenotypes, and therefore the overall population stratification will not confound the association test as it does for conventional case-control tests.

A total of 153,348 bi-allelic SNPs and 15,468 InDels were detected across the core genome for the training cohort GAS strains. Among them, 138,907 (82.3%) fall into coding regions and 63,321 (37.5%) cause translational changes (non-synonymous SNPs and frame-shift InDels) (Figure S3a). The gene-wise profiling of the mutations reveals a quasi-normal distribution across the core genes with an average number of 80 mutations per gene, corresponding to an average mutation ratio of 8.6% per gene (Figure S3b). The majority of the mutations have a minor allele frequency (MAF) < 10% in all phenotypes, but the distribution pattern of MAF varies among the phenotypes (Figure S3c). Further inspection of the pairwise correlation of MAF between phenotypes shows that DSTI is significantly different from other phenotypes, and therefore, was removed from the association test (Figure S3d-S3i). We used the remaining seven phenotypes (INV, PHY, SSTI, RF, SF, SP and PM) for downstream analysis. We then took a strategy of case-case association analysis to identify genetic signatures differentiating pair-wise phenotypes of GAS by constructing ten binary comparison groups for association tests, i.e., INV-PHY, INV-SSTI, INV-PHYSSTI, INV-RF, PHY-SSTI, PHY-RF, RF-SSTI, PM-SF, PM-SP, and SF-SP.

New InDel-based and CDS-based GWAS methods enable detection of missing association between GAS phenotypes and genotypes

In order to perform association tests for a broader range of genetic variants in addition to SNPs, we developed new methods IndelWinScan for InDel-based GWAS and CDS2GWAS for CDS-based GWAS (Fig. 2a), which allowed association tests for InDel mutational status of nonCDS fragments and the coding status of CDS, respectively. Our methods treated nonCDS or CDS as genetic units and evaluated the overall genetic status of the units caused by integrative effects of SNPs and InDels within the units, and then transformed the genetic status to categorical genotypes, i.e., intact vs. inactivated or intact vs. mutated (see Methods & Materials). The transformed genotypes were finally used for testing associations with the phenotypes of interest. For those genes without being inactivated, we employed the traditional SNP-based GWAS analysis to identify the associations between SNPs in those genes and phenotypes of interest.

Using these newly developed methods, we performed GWAS analysis in three modalities across the whole genomes of GAS, namely, SNP-based GWAS, InDel-based GWAS, and CDS-based GWAS for each of the ten paired groups of phenotypes we constructed. The association results for the ten groups are summarized in Fig. 2b. It is seen that the highest number of genotype-phenotype associations was detected from the coding SNPs in all tested groups, indicating that SNPs in the coding regions are still the major source of genetic variations. The InDel mutational fragments in nonCDS regions are the second leading source of genetic variations with phenotypic associations, highlighting the importance of investigating mutations in these regions. The truncated CDS genes account for the lowest number of associations.

In addition to the genomic region-specific distribution of the variants, a clear phenotype-dependent pattern was also observed. For instance, the INV-containing comparison groups (INV-PHY, INV-SSTI and INV-PHYSSTI) yield the highest number of truncated CDS significantly associated with phenotype differentiation, suggesting frequent alterations of gene coding status in invasive GAS strains in comparison with non-invasive strains.

To further assess the loci associated with the phenotypes, we took the comparison groups INV-SSTI as an example to visualize the variants across the genome with the Manhattan plot (Fig. 3, the full set of Manhattan plots can be found in Figure S4) and then further validated the most significant candidates in the validation cohort of GAS strains (Fig. 4). The most significant variants from the SNP-based GWAS were detected in the gene nga, encoding NAD glycohydrolase. Seven of the variants were confirmed between INV-SSTI in the validation cohort, including two non-coding SNPs in the promoter spacer region (T-22G at 150429 and C-18T at 150433) and four non-synonymous SNPs in the coding region (Ala99Val, Arg103His, Gly136Arg, Met195Ile, and Leu199Ile) with high association significance (P < 10^− 3, Fig. 4a). The two SNPs in the nga promoter at -22 and − 18 have been reported to be related with enhanced virulence in animal models of pharyngitis and necrotizing fasciitis [32].

For InDel-based GWAS, 244 mutated non-coding fragments showed significant association in the testing for the group INV-SSTI. Those fragments are located at the 5’ upstream regions of 37 genes (with an average distance of 198 bp) or the 3’ downstream regions of 32 genes (with an average distance of 214 bp) (Fig. 3b and Figure S5). Among these fragments, the two most significant fragments associated with INV are in the promoter region of the small subunit ribosomal RNA (SSU rRNA) located at around 1.59 Mb (AP53_1607) and 16 kb (AP53_15) of the reference genome, respectively (P < 10^-4) (Fig. 4b and 4c). The two fragments were further verified in the validation cohort, which contain multiple InDels, including one in the UP element of the SSU rRNA promoter region [33]. Previous studies showed that the UP element can enhance rRNA transcription upon bound by the RNA polymerase α subunit in both Escherichia coli and Bacillus subtilis, although to different extents [34, 35].

Notably, for the CDS-based GWAS, two known virulence regulator gene, covS and rocA, were identified to be among the most significant genes whose inactivated coding status is associated with INV relative to SSTI (Fig. 3c). Although the inactivation genotype of covS has been experimentally determined to be associated with hypervirulence and invasiveness in several GAS isolates [13–16], our study first identifies the association on the population level by leveraging the advantage of our new CDS-based GWAS method, which allows testing the gene coding status regardless of the specific variant forms in the gene. This association was confirmed in the validation cohort (Fig. 4d-e). Specifically, we detected 12 different forms of allelic variants, including SNPs, insertions, and deletions in the covS coding region independently occurred in 26 GAS strains (Fig. 4d). The 12 variants all led to covS inactivation and when tested collectively as a gene-inactivation genotype, showed highly significant association with GAS invasiveness (β = 2.3, P = 0.002). Similarly, we detected three different forms of inactivating mutations in the gene rocA from seven GAS strains. By treating them collectively as a gene-inactivation genotype, we identified the significant association between the inactivated rocA with GAS invasiveness (β = 2.1, P = 0.017) (Fig. 4e). Previous experiments demonstrated that truncated rocA or suppressed expression of rocA resulted in enhanced GAS survival in human blood in M18, M3, and M89 GAS isolates [36–38]. The successful identification of the association of covS/rocA inactivation with GAS invasiveness demonstrates the capability of our new method of GWAS in detecting key associations which have been missed by other methods.

The variants in different genomic partitions are involved in distinct biological functions

In order to characterize the pattern of genotype-phenotype associations at different genomic regions, we made genomic partitioning for the associated variations detected in the three GWAS modalities into four categories based on their genic locations and functional impacts, i.e., the 5’ upstream non-coding regions of genes (5’-nonCDS), the 3’ downstream non-coding regions of genes (3’-nonCDS), intact genes without being truncated (intactCDS), truncated CDS by any kind of mutations (trunCDS) (Fig. 5a). To reduce the influence of redundant signals, we only selected the genetically independent and functionally effective variants using a pruning and thresholding strategy based on linkage disequilibrium (LD) and functional annotation [39] (Figure S6, See Methods & Materials). As a result, we obtained 1,361 variants linked with 783 genes across the four partitions (Additional File 5). On the one hand, 251 (32.0%) genes are linked with variants from more than one partition and 496 (63.3%) have variants significant in more than one binary-phenotype tests (Fig. 5b). On the other hand, all 10 binary-phenotype tests detected associated variants from at least two partitions and at least two genes (Figure S7). The genotype-phenotype association relationships clearly reveal a pleiotropic and polygenic pattern.

To obtain functional insights of the variants significant in the association tests in distinct genomic partitions, we classified the 783 genes linking with associated variants into 30 functional classes (Fig. 5c and Additional File 6). The distribution of the genes in the functional classes is highly dependent on the genomic partitions of the variants and the binary-phenotype tested. Specifically, the variants significant in the phenotype pairs INV-SSTI, INV-PHYSSTI, RF-SSTI, PM-SF, and SF-SP are more frequently localized in the genes involved in metabolism, such as “Carbohydrates”, “DNA Metabolism”, “Protein Metabolism”, and “RNA Metabolism” (P < 0.05), whereas the variants significant in the phenotype pairs INV-RF, PHY-RF, and PM-SP are sporadic in the metabolic functional classes, suggesting higher metabolic similarities between the three phenotype pairs. The associated variants in the partition of 5’-nonCDS are more significantly enriched in the genes involved in “Virulence, Disease and Defense”, and the associated SNPs in the partition of intactCDS are enriched in several categories, such as “Virulence, Disease and Defense”, “Protection System”, “Iron Acquisition and Metabolism”, and “Carbohydrates”, indicating multi-faceted influences of the coding SNPs on a wide range of functions in GAS. The translational truncation of coding genes (trunCDS) avoids most of the functional classes but is predominantly enriched in the category “Regulation and cell signaling”. The truncated genes in this category include not only the well-known virulence regulator covS, but also several other less-studied regulators, such as rocA, rivR, mutR, and luxR. It indicates that the transcriptional shift caused by transcriptional regulators is a common mechanism for GAS to switch the pathogenicity or phenotypes [40, 41].

Genetic architecture of GAS phenotypes is a hybrid model

The polygenic nature of GAS phenotypes prompted us to investigate the detailed genetic architecture of the phenotypes. By examining the relationship between MAF and effect size of the genetic variants detected in GWAS, our data supports a multifaceted genetic architecture of GAS phenotypes, which is manifested as a combination of a long tail of common variants (MAF ≥ 5%) with low effect sizes and a small proportion of low-frequency variants (1% ≤ MAF < 5%) with large effect sizes (Fig. 6a-d). Of particular interest is the partition of trunCDS, where the genetic architecture is dominated by low-frequency variants of large effect size, signifying a trend of lower population frequency and elevated impact of gene inactivation on phenotypic differentiation of GAS relative to the non-inactivating variants. A notable example in this partition is rivR, a regulator of virulence factor expression, which has an effect size as high as e^3.3~27.1 and MAF ~ 2.1%. Interestingly, the covS inactivation in this partition also has a high effect size of e^3.0~20.1, but the MAF ~ 7.8% at a modest level. It provides the evidence for the presence of high-effect variants segregating at an appreciable frequency in polygenic phenotypes of GAS.

Overall, the genetic architecture of the complex phenotypes of GAS is characterized by highly polygenic variants covering a wide spectrum of allelic frequency and effect size along with a few close-to-common variants with large effect sizes.

To further delineate the functional implications of the polygenic genetic architecture of GAS, we evaluated the association on the gene level and function level. We first constructed the bipartite network graph between the eight phenotypes and the genes in the four genomic partitions (Fig. 6e-h). A phenotype and a gene are connected by a positive/negative edge if that gene is linked with variants positively/negatively associated with that phenotype. Each of the eight phenotypes is associated with multiple genes across the four partitions and the gene number varies in a wide range from 6 to 148 depending on the specific phenotypes (Fig. 6i and Additional File 7). Regardless of the genomic partitions, PHY has the most positively associated genes and SSTI has the most negatively associated genes, indicating substantial genetic differences between these two major symptomatic phenotypes. In the partition of trunCDS, the phenotype INV is associated with multiple inactivated genes in both positive and negative direction with covS the most significant one (P = 5.0 x 10^− 5).

By integrating the contributions of multiple genes involved in the same functional category, we then constructed the bipartite network between the eight phenotypes and the functional categories (Figure S8). The overall network reflects the complex multifaceted etiology of GAS phenotypes and pathogenesis. The phenotypes are specifically associated with certain functional categories in each partition. For instance, the phenotype INV is most significantly associated with the functional category “Virulence, Disease and Defense” via the SNPs in the partition of intactCDS (P = 5.1 x 10^− 16, weight = 15.3) and via the variants in the partition of 5’-nonCDS (P = 4.0 x 10^− 18, weight = 17.4). Its association with the category “Regulation and Cell signaling” is prominent in the partition of trunCDS (P = 2.2 x 10^− 7, weight = 6.65). The strong association between INV and these two functional categories highlights the key functional modifications and the way the functions were modified underlying GAS invasiveness. For PHY and SSTI which have the highest number of positive and negative gene-level associations, respectively, the two phenotypes exhibit highly significant and opposite associations with several common functional categories (P < 10^− 8), such as “Carbohydrates” (P < 2.5 x 10^− 11), “DNA Metabolism” (P < 4.0 x 10^− 9), “RNA Metabolism” (P < 3.1 x 10^− 9), “Protein Metabolism” (P < 10^− 8), “Virulence, Disease and Defense” (P < 1.6 x 10^− 12), “Cofactors, Vitamins, and Prosthetic Groups” (P < 1.3 x 10^− 5) in the partition of intactCDS or 5’-nonCDS. Given that GAS strains of PHY and SSTI phenotype have been known to diverge in two virulence-related gene clusters, i.e., mga regulon and pilus assembly, our data here provides complementary genetic factors contributing to the differences between these two highly prevalent GAS phenotypes. The data supports a paradigm of phenotypic differentiation of GAS determined by a combination of high divergence in several key genes and many small-effect variants in genes responsible for a wide spectrum of cellular functions, which together shape final pathogenesis, survival and adaptation capacity of GAS.

Pleiotropicity of the associated genes

As a reciprocal side of polygenicity, we next investigated the pleiotropicity of the associated genes across the eight phenotypes. It is manifested by an average of ~ 38.6% of the genes associated with at least two phenotypes and 5.3% associated with at least four phenotypes (Fig. 7a-d and Additional File 7). Further functional enrichment analysis showed that the genes with pleiotropicity ≥ 4 are more significantly enriched in the category “Virulence, Disease and Defense” (in the genomic partition of intactCDS and 5’-nonCDS) and “Iron Acquisition and Metabolism” (in the partition of intactCDS), in clear contract to the functional enrichment pattern for the genes with pleiotropicity < 4, which expands a wide spectrum of categories, such as “Regulation and Cell signaling”, “Iron Acquisition and Metabolism”, and “Cofactors, Vitamins, and Prosthetic Groups” (Fig. 7e) [42, 43]. The distinct functional enrichment pattern of polygenicity groups highlights the critical roles of the genes in the two functional categories “Virulence, Disease and Defense” and “Iron Acquisition and Metabolism” in diversifying GAS phenotypes or pathogenesis.

The pleiotropicity of the associated genes is further illuminated by 26 genes harboring distinct non-synonymous SNPs within the same genes but associated oppositely with the same phenotype (Additional File 7). Among them, the most notable gene is nga which is associated with all eight phenotypes via ten distinct variants from three genomic partitions, i.e., the 5’-nonCDS, 3’-nonCDS, and intactCDS (Fig. 7f and Additional File 8).

The genetic variants from three GWAS modalities explained a large proportion of GAS phenotypic variance

Next, we assessed the phenotypic variance explained by the polygenic variants detected in each phenotype-pair test. That will help to estimate the potential of the variants for GAS phenotype differentiation or disease-causal prediction. In order to minimize the influence of co-linear variants, we first selected the variants without mutual co-linearity within partitions using the LASSO regression method (Figure S9 and Additional File 9), and then calculated the contribution of the variants to the phenotypic variance using the coefficients of determinant measure (R²) on the liability scale [31, 44]. The relative proportions of variance explained by the variants from different GWAS modalities are highly dependent on the phenotypes tested (Fig. 8a). For the phenotype pairs involving INV, i.e., INV-PHY, INV-PHYSSTI, INV-SSTI, and INV-RF, the variants from the SNP-based GWAS (SNP-coding + SNP-nonCoding) explained 34%-48% of the pair-wise phenotypic variance and the gene inactivations from the CDS-based GWAS explained an additional 10% (4%-15%) of the variance, supporting the roles of gene inactivation for GAS invasiveness transition. Again for the four phenotype pairs, the variants from the InDel-based GWAS explained another 0%-4% of the variance. The combination of the gene inactivation from CDS-based GWAS and InDels from InDel-based GWAS together contributed additional 7%-16% of the phenotypic variance, pushing the total average explained variance above 50% (47%-57%). The results show that our new method of GWAS is able to significantly narrow the gap between GAS phenotypes and genotypes.

An exception is the phenotype pair PHY-RF, where all variants only explained a total of 27% variance, implying that many other factors may contribute to the differences between PHY and RF, such as epigenetic modifications and host interactions. For the other phenotype pairs PM-SF, SF-SP, and PM-SP, the explained variance did not increase by additional types of variants. This scenario is partially due to the small sample size of the three phenotypes which biased the estimation of effect sizes and association significance.

Development of a GWAS polygenic score model for phenotypic differentiation of GAS

The large proportion of explained variance by the variants from multiple GWAS modalities motivated us to explore the capability of the variants as genetic markers for GAS phenotype stratification or disease prediction. The polygenic nature of the phenotypes enabled us to develop a GWAS polygenic score (GPS) as a surrogate of multiple genetic markers. GPS is a composite measure of genetic mutation burden over polygenic loci for single individuals and has been frequently used for human disease prediction or risk stratification [45–47]. Here, by leveraging the polygenic nature of GAS phenotypes, we developed GPS-based prediction models for GAS phenotype differentiation by summing over the independent variants from three GWAS modalities (Figure S9 and Additional File 9 (see Methods & Materials). We then evaluated the performance of the GPS models in discriminating GAS phenotypes in the validation cohort of 710 GAS strains using the Receiver Operating Characteristic (ROC) analysis (Additional File 10). The strains of the validation cohort have the similar emm serotype distribution and time coverage to those of the training cohort (Fig. 8b and 8c). The ROC analysis showed that the GPS calculated from each of the three GWAS modalities alone (i.e., SNP-based, CDS-based, and InDel-based GWAS) generated models with an average AUC of 77.5%, 64.4%, and 69.4%, respectively (Fig. 8d). Noticeably, incorporating the CDS-based GWAS made AUC increase by 3.8% and combining all three GWAS modalities made AUC increase by 4.4%, leading to a high discriminatory level above 80% (Fig. 8d). Our results demonstrate the considerable improvement in the performance of the GPS prediction models by adding CDS-based and InDel-based GWAS markers. This high level of discriminatory ability of the integrated GPS score facilitates its promising role in phenotype stratification for this highly versatile organism.

A low-end exception in the AUC values occurred for the phenotype pair PHY-RF, which has an AUC of 74.4%. This modest value of AUC for PHY-RF agrees well with the small proportion of explained variance (27%) by all variants between the two phenotypes, further supporting the possibility of other contributing factors than genetic mutations to the phenotypic difference between PHY-RF.

GAS is a very unique pathogen. It has a genome of length shorter than 2 Mbp, but is capable of causing a variety of human diseases. Its disease causal capacity is largely attributed to the highly recombining genome, high level of genetic variation, numerous virulence factors, and intricate regulatory circuits. Over the decades, numerous genetic factors as well as variations in these genetic factors have been identified to be key determinants for the pathogenesis or disease phenotypes. The virulence genes are the most notable group of factors associated with the pathogenesis, such as M protein, cysteine protease (speB), and streptolysin O (slo), streptolysin S (sagA), IL8 peptidase (scpC), streptokinase (ska), C5a peptidase (scpA), esterase (sse), hyaluronic acid capsule (hasABC), and iron acquisition [2, 5, 6]. The two-component regulatory system covRS is the most prominent regulatory gene determining GAS invasiveness. Later, many metabolic genes were also found to contribute to the virulence of GAS, such as ccpA [48, 49] and MtsR [50, 51]. However, these genes were usually studied individually in a specific strain or a specific type of diseased mouse model. They have not been systematically investigated in a population level in such a manner as to clarify the many-to-many relationships between genetic factors and phenotypes or reveal the genetic underpinnings of the phenotypic versatility of GAS. Therefore, the early diagnosis of GAS has been difficult. Our report is the first such study to delineate these issues by leveraging our newly developed GWAS tool for identifying the genetic architecture and full spectrum of genetic variants including InDels and gene inactivation in a bacterial species.

Our data clearly revealed a high level of polygenicity of the GAS phenotypes and pleiotropicity of the causal genes, whereby each disease phenotype is associated with multiple causal mutations/genes and the single gene may contain distinct mutations associated to different phenotypes. These two properties form the foundation of complex relationship between numerous causal genes and versatile phenotypes of GAS. The most polygenic phenotype is the skin-related SSTI, which is associated with more than hundreds of genes spanning a variety of functions, including metabolism and virulence. For the least polygenic phenotype INV, although covS inactivation is the most significant causal variation, this phenotype is still manifested to be associated with multiple lines of gene changes from different genomic partitions, particularly the SNPs from 87 coding genes involved in diverse functions, such as virulence, carbohydrate metabolism, protein metabolism, nucleotide metabolism among others. Although the genetic engineering experiments have repeatedly demonstrated that mutation-caused covS/rocA inactivation alone is sufficient to derive an isogenic invasive GAS phenotype [6, 13, 37, 52], our data suggests that naturally evolved invasive strains may have undergone multiple mutations in a series of genes in addition to covS, the optimized combination of which ultimately determined specific disease phenotypes and contributed to adaptation capabilities.

However, we observed that the manner in which the mutations influence the genes and functions is heterogeneous across genomic partitions. For instance, the genetic variants in the 5’ non-coding regions are predisposed to influence the virulence-related functions, and the coding SNPs within genes tend to make changes in many different functions, including virulence, carbohydrate metabolism, and iron acquisition, whereas the truncating mutations avoid most of the functional classes but predominantly inactivate genes responsible for regulation. This observation echoes the current consensus on the regulator genes like covS or rocA, that the solely inactivation of either one can modulate the expression of a large proportion of genes in GAS and consequently result in invasive transition [53]. The covS/rocA inactivations are prominent examples of cost-effective and large-effect variations occurred in GAS, which are key players in determining invasiveness. It is interesting that inactivation of regulator genes is not the only variant type that has large effects on phenotypes. The pharyngitis (PHY), the most common symptomatic phenotype of GAS, is associated with a SNP in the 5’ non-coding region of nga (roughly 300bp upstream of the nga promoter) with a large effect size (e^3.0~36.6), in addition to many other variants across all genomic partitions with low effect sizes. It is reasonable to speculate that the two gene clusters, i.e., mga regulon and pilus assembly, known to highly diverge and differentiate between PHY and other types [54, 55], should also have large effect on phenotype differentiation. Taken together, our data revealed a genetic architecture of GAS phenotypes as a combination of many low-effect common variants and a small proportion of high-effect low-frequency variants as demonstrated in the curve of effect size versus MAF. We noticed an appreciable level of MAF for covS inactivation (~ 7.8%) in our studied dataset probably induced from frequent genetic engineering of covS-inactivated mutant strains or intentional isolation of covS-inactivated pathogenic strains for research purpose. It is expected that the true frequency of covS inactivation in natural population should be at the level of rare variants.

Consistent with the large effect size of gene inactivations, our calculation showed that the gene inactivations (or the genetic variants in the partition trunCDS) only occupy 8.7% of the total detected variants, however, they are able to explain 31%-56% of the phenotypic variance for the phenotype pairs INV-PHY, INV-PHYSSTI, and INV-SSTI, comparable to or higher than the proportion explained by the coding SNPs, which comprise 30.2% of the total variants. When combining the contributions of all variant types, the average explained variance reached above 50%. Although there has been no such explained variance reported in GAS and other bacterial organism, this proportion is comparable to that for the human complex trait height (45%-52%) and higher than that for the systolic blood pressure (22%-25%) [56–58]. It indicates that our newly developed GWAS methods have been able to detect genetic variants associated with GAS phenotypes to a relatively full extent. We noted that the total proportion of explained variance did not proportionally increased by combining multiple types of variants, specifically the gene inactivations explaining an additional 10% (4%-15%) of the variance and the InDels explaining another 0%-4% of the variance, significantly lower than that explained by each alone. Possible explanations include the existence of gene-gene interactions, gene-host interactions, and epigenetic modifications. Future studies of the influences of these factors to phenotypic variance are warranted.

Based on the current level of explained variance, the genetic variants allowed us to develop the composite GPS models for phenotype prediction or risk stratification. Incorporating the gene inactivations into the models improved the AUC by 3.8% and combining all types of variants led to the AUC above 80%, demonstrating the potential of the GPS serving as biomarkers for guiding early diagnosis and curtailing progression of GAS diseases.

Even with the promising results, we realize that there are several issues to be addressed here and improved upon in our future study. First, the structural variations were not incorporated into our GWAS method due to the sparsity of this variant type in GAS genomes. Secondly, gene-gene and gene-host interactions were not considered in the association tests. We expect that new methods for modeling the influence of these interactions will provide novel insights of mechanisms underlying GAS phenotypic complexity. Thirdly, the LD between variants from different GWAS modalities was not used for variant pruning. It may also contribute to the large shared proportion of explained phenotypic variance between variants from SNP-based GWAS and CDS-based GWAS.

We have provided a new GWAS tool and analysis framework for comprehensively identifying genetic variants associated with pathogenic phenotypes. It allowed for revealing the complete genetic architecture and developing risk stratification models in the highly versatile pathogen of GAS. The new tool can be broadly applied in other bacterial species of complex traits or disease phenotypes.

Competing Interest

The authors declare no conflicts of interest.

Ethics approval and consent to participate

No ethical approval was required as no human or animal experiments were involved.

Consent for publication

The authors declare that there are no potential competing interests.

Availability of data and materials

The data used in this study is provided as additional files. The package and source code of compGWAS is available at: https://github.com/BaoCodeLab/compGWAS under the MIT license.

Acknowledgements

Not applicable.

Funding

The work was supported by a faculty Initiation Grant Program.

Authors’ contributions

YJB and FJC conceived and supervised the study. PYW and YJB performed the analysis and developed the software. ZL performed experimental studies. PYW, ZL, and ZSC performed data mining and curation. YB and FJC wrote the manuscript. All authors revised and approved the final manuscript.

Carapetis JR, Steer AC, Mulholland EK, Weber M. The global burden of group A streptococcal diseases. Lancet Infect Dis. 2005;5(11):685-94.
Walker MJ, Barnett TC, McArthur JD, Cole JN, Gillen CM, Henningham A, et al. Disease manifestations and pathogenic mechanisms of Group A Streptococcus. Clin Microbiol Rev. 2014;27(2):264-301.
Logan LK, McAuley JB, Shulman ST. Macrolide treatment failure in streptococcal pharyngitis resulting in acute rheumatic fever. Pediatrics. 2012;129(3):e798-802.
Misiakos EP, Bagias G, Patapis P, Sotiropoulos D, Kanavidis P, Machairas A. Current concepts in the management of necrotizing fasciitis. Frontiers in surgery. 2014;1:36.
Barnett T, Indraratna A, Sanderson-Smith M. Secreted Virulence Factors of Streptococcus pyogenes. In: Ferretti JJ, Stevens DL, Fischetti VA, editors. Streptococcus pyogenes: Basic Biology to Clinical Manifestations. Oklahoma City (OK): University of Oklahoma Health Sciences Center © The University of Oklahoma Health Sciences Center.; 2022.
Cole JN, Barnett TC, Nizet V, Walker MJ. Molecular insight into invasive group A streptococcal disease. Nat Rev Microbiol. 2011;9(10):724-36.
Power RA, Parkhill J, de Oliveira T. Microbial genome-wide association studies: lessons from human GWAS. Nat Rev Genet. 2017;18(1):41-50.
Kachroo P, Eraso JM, Beres SB, Olsen RJ, Zhu L, Nasser W, et al. Integrated analysis of population genomics, transcriptomics and virulence provides novel insights into Streptococcus pyogenes pathogenesis. Nat Genet. 2019;51(3):548-59.
Bao Y-J, Shapiro BJ, Lee SW, Ploplis VA, Castellino FJ. Phenotypic differentiation of Streptococcus pyogenes populations is induced by recombination-driven gene-specific sweeps. Sci Rep. 2016;6:36644.
Davies MR, McIntyre L, Mutreja A, Lacey JA, Lees JA, Towers RJ, et al. Atlas of group A streptococcal vaccine candidates compiled using large-scale comparative genomics. Nat Genet. 2019;51(6):1035-43.
Mosquera-Rendón J, Moreno-Herrera CX, Robledo J, Hurtado-Páez U. Genome-Wide Association Studies (GWAS) Approaches for the Detection of Genetic Variants Associated with Antibiotic Resistance: A Systematic Review. Microorganisms. 2023;11(12):2866.
San JE, Baichoo S, Kanzi A, Moosa Y, Lessells R, Fonseca V, et al. Current Affairs of Microbial Genome-Wide Association Studies: Approaches, Bottlenecks and Analytical Pitfalls. Front Microbiol. 2019;10:3119.
Sumby P, Whitney AR, Graviss EA, DeLeo FR, Musser JM. Genome-Wide Analysis of Group A Streptococci Reveals a Mutation That Modulates Global Phenotype and Disease Specificity. PLoS Pathog. 2006;2(1):e5.
Mayfield JA, Liang Z, Agrahari G, Lee SW, Donahue DL, Ploplis VA, et al. Mutations in the control of virulence sensor gene from Streptococcus pyogenes after infection in mice lead to clonal bacterial variants with altered gene regulatory activity and virulence. PLoS One. 2014:in press.
Bao Y, Liang Z, Booyjzsen C, Mayfield JA, Li Y, Lee SW, et al. Unique genomic arrangements in an invasive serotype M23 strain of Streptococcus pyogenes identify genes that induce hypervirulence. J Bacteriol. 2014;196(23):4089-102.
Garcia AF, Abe LM, Erdem G, Cortez CL, Kurahara D, Yamaga K. An insert in the covS gene distinguishes a pharyngeal and a blood isolate of Streptococcus pyogenes found in the same individual. Microbiology. 2010;156(Pt 10):3085-95.
Lees JA, Galardini M, Bentley SD, Weiser JN, Corander J. pyseer: a comprehensive tool for microbial pangenome-wide association studies. Bioinformatics. 2018;34(24):4310-2.
Lees JA, Vehkala M, Välimäki N, Harris SR, Chewapreecha C, Croucher NJ, et al. Sequence element enrichment analysis to determine the genetic basis of bacterial phenotypes. Nat Commun. 2016;7:12797.
Aun E, Brauer A, Kisand V, Tenson T, Remm M. A k-mer-based method for the identification of phenotype-associated genomic biomarkers and predicting phenotypes of sequenced bacteria. PLoS Comput Biol. 2018;14(10):e1006434.
Earle SG, Wu CH, Charlesworth J, Stoesser N, Gordon NC, Walker TM, et al. Identifying lineage effects when controlling for population structure improves power in bacterial association studies. Nature microbiology. 2016;1:16041.
Bao Y-J, Liang Z, Mayfield JA, Donahue DL, Carothers KE, Lee SW, et al. Genomic characterization of a pattern D Streptococcus pyogenes emm53 isolate reveals a genetic rationale for invasive skin tropicity. J Bacteriol. 2016;198:1712-24.
Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22(3):568-76.
Letunic I, Bork P. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 2021;49(W1):W293-w6.
Stirling JR, Zakynthinaki M. The Point of Maximum Curvature as a Marker for Physiological Time Series. Journal of Nonlinear Mathematical Physics. 2008;15(3):396-406.
Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005;21(2):263-5.
Zhang S, Chen H-S, Pfeiffer RM. A combined p-value test for multiple hypothesis testing. Journal of Statistical Planning and Inference. 2013;143(4):764-70.
Heard NA, Rubin-Delanchy P. Choosing between methods of combining $p$-values. Biometrika. 2018;105(1):239-46.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498-504.
McKelvey RD, Zavoina W. A statistical model for the analysis of ordinal level dependent variables. The Journal of Mathematical Sociology. 1975;4(1):103-20.
Heinze G, Schemper M. A solution to the problem of separation in logistic regression. Statistics in medicine. 2002;21(16):2409-19.
Lee SH, Goddard ME, Wray NR, Visscher PM. A better coefficient of determination for genetic profile analysis. Genetic epidemiology. 2012;36(3):214-24.
Zhu L, Olsen RJ, Nasser W, Beres SB, Vuopio J, Kristinsson KG, et al. A molecular trigger for intercontinental epidemics of group A Streptococcus. J Clin Invest. 2015;125(9):3545-59.
Shin Y, Qayyum MZ, Pupov D, Esyunina D, Kulbachinskiy A, Murakami KS. Structural basis of ribosomal RNA transcription regulation. Nat Commun. 2021;12(1):528.
Krásný L, Gourse RL. An alternative strategy for bacterial ribosome synthesis: Bacillus subtilis rRNA transcription regulation. Embo j. 2004;23(22):4473-83.
Maeda M, Shimada T, Ishihama A. Strength and Regulation of Seven rRNA Promoters in Escherichia coli. PLoS One. 2015;10(12):e0144697.
Lynskey NN, Goulding D, Gierula M, Turner CE, Dougan G, Edwards RJ, et al. RocA truncation underpins hyper-encapsulation, carriage longevity and transmissibility of serotype M18 group A streptococci. PLoS Pathog. 2013;9(12):e1003842.
Miller EW, Danger JL, Ramalinga AB, Horstmann N, Shelburne SA, Sumby P. Regulatory rewiring confers serotype-specific hyper-virulence in the human pathogen group A Streptococcus. Molecular microbiology. 2015;98(3):473-89.
Zhu L, Olsen RJ, Horstmann N, Shelburne SA, Fan J, Hu Y, et al. Intergenic Variable-Number Tandem-Repeat Polymorphism Upstream of rocA Alters Toxin Production and Enhances Virulence in Streptococcus pyogenes. Infect Immun. 2016;84(7):2086-93.
Khera AV, Chaffin M, Aragam KG, Haas EM, Roselli C, Choi SH, et al. Genome-wide polygenic scores for common diseases identify individuals with risk equivalent to monogenic mutations. Nat Genet. 2018;50(9):1219-24.
Lee LW, Mapp AK. Transcriptional switches: chemical approaches to gene regulation. J Biol Chem. 2010;285(15):11033-8.
Beres S, Kachroo P, Nasser W, Olsen R, Zhu L, Flores A, et al. Transcriptome remodeling contributes to epidemic disease caused by the human pathogen Streptococcus pyogenes. mBio. 2016;7(3):e00403-16.
Ge R, Sun X. Iron acquisition and regulation systems in Streptococcus species. Metallomics : integrated biometal science. 2014;6(5):996-1003.
Fisher M, Huang YS, Li X, McIver KS, Toukoki C, Eichenbaum Z. Shr is a broad-spectrum surface receptor that contributes to adherence and virulence in group A streptococcus. Infect Immun. 2008;76(11):5006-15.
Zhu H, Zhou X. Statistical methods for SNP heritability estimation and partition: A review. Computational and structural biotechnology journal. 2020;18:1557-68.
Torkamani A, Wineinger NE, Topol EJ. The personal and clinical utility of polygenic risk scores. Nat Rev Genet. 2018;19(9):581-90.
Lennon NJ, Kottyan LC, Kachulis C, Abul-Husn NS, Arias J, Belbin G, et al. Selection, optimization and validation of ten chronic disease polygenic risk scores for clinical implementation in diverse US populations. Nat Med. 2024;30(2):480-7.
Lewis ACF, Perez EF, Prince AER, Flaxman HR, Gomez L, Brockman DG, et al. Patient and provider perspectives on polygenic risk scores: implications for clinical reporting and utilization. Genome Med. 2022;14(1):114.
Paluscio E, Watson ME, Jr., Caparon M. CcpA Coordinates Growth/Damage Balance for Streptococcus pyogenes Pathogenesis. Sci Rep. 2018;8(1):14254.
DebRoy S, Aliaga-Tobar V, Galvez G, Arora S, Liang X, Horstmann N, et al. Genome-wide analysis of in vivo CcpA binding with and without its key co-factor HPr in the major human pathogen group A Streptococcus. Molecular microbiology. 2021;115(6):1207-28.
Merchant AT, Spatafora GA. A role for the DtxR family of metalloregulators in gram-positive pathogenesis. Molecular oral microbiology. 2014;29(1):1-10.
Do H, Makthal N, VanderWal AR, Saavedra MO, Olsen RJ, Musser JM, et al. Environmental pH and peptide signaling control virulence of Streptococcus pyogenes via a quorum-sensing pathway. Nat Commun. 2019;10(1):2586.
Bernard PE, Kachroo P, Zhu L, Beres SB, Eraso JM, Kajani Z, et al. RocA Has Serotype-Specific Gene Regulatory and Pathogenesis Activities in Serotype M28 Group A Streptococcus. Infect Immun. 2018;86(11):e00467-18.
Vega LA, Malke H, McIver KS. Virulence-Related Transcriptional Regulators of Streptococcus pyogenes. In: Ferretti JJ, Stevens DL, Fischetti VA, editors. Streptococcus pyogenes: Basic Biology to Clinical Manifestations. Oklahoma City (OK): University of Oklahoma Health Sciences Center; 2022.
Bessen DE, Lizano S. Tissue tropisms in group A streptococcal infections. Future microbiology. 2010;5(4):623-38.
Kreikemeyer B, Gámez G, Margarit I, Giard J-C, Hammerschmidt S, Hartke A, et al. Genomic organization, structure, regulation and pathogenic role of pilus constituents in major pathogenic Streptococci and Enterococci. Int J Med Microbiol. 2011;301(3):240-51.
Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42(7):565-9.
Akiyama M, Ishigaki K, Sakaue S, Momozawa Y, Horikoshi M, Hirata MM, K.Ikegawa, S.Takahashi, A., et al. Characterizing rare and low-frequency height-associated variants in the Japanese population. Nat Commun. 2019;10(1):4393.
Aslibekyan S, Wiener HW, Wu G, Zhi D, Shrestha S, de Los Campos G, et al. Estimating proportions of explained variance: a comparison of whole genome subsets. BMC proceedings. 2014;8(Suppl 1 Genetic Analysis Workshop 18Vanessa Olmo):S102.

No competing interests reported.

AdditionalFile1.xlsx
Additional File 1.xlsx: Catalog of the 2,598 GAS strains and clinical information collected for this study.
AdditionalFile2.xlsx
Additional File 2.xlsx: Catalog of the 165 GAS strains used for constructing the core genome.
AdditionalFile3.xlsx
Additional File 3.xlsx: List of 104 genes excluded for regression modeling in GWAS analysis.
AdditionalFile4.xlsx
Additional File 4.xlsx:Catalog of the 790 GAS strains used as the training cohort.
AdditionalFile5.xlsx
Additional File 5.xlsx:List of 1,361 genetic variants associated with GAS phenotypes in the four genomic partitions.
AdditionalFile6.xlsx
Additional File 6.xlsx:Frequency distribution of the 783 genes linking with 1361 variants among the functional categories.
AdditionalFile7.xlsx
Additional File 7.xlsx:The gene-phenotype matrix presentation of the polygenicity of GAS phenotypes and pleoitropicity of GAS genes in the four genomic partitions. The association significance in –log10(p-value) was shown. The 26 genes harboring distinct non-synonymous SNPs within the same genes but associated oppositely with the same phenotype are also listed.
AdditionalFile8.xlsx
Additional File 8.xlsx:The detailed information of the ten genetic variants in the gene nga, including GWAS summary statistics and allele frequencies.
AdditionalFile9.xlsx
Additional File 9.xlsx:The list of selected genetic variants of GAS used for evaluation of explained variance and GPS.
AdditionalFile10.xlsx
Additional File 10.xlsx:Catalog of the 710 GAS strains used as the testing cohort.
FigureS1.pdf
Supplementary Figures
Figure S1: Heatmap presentation of the mapping coverage of the 1584 core genes among the 165 complete GAS genomes. The emm type of the genomes is shown at the top. The genes were hierarchically clustered.
FigureS2.pdf
Figure S2: The phylogenetic tree of the training cohort of 790 GAS strains built with the core genomic sites. Annotations are shown from the innermost to the outermost circles: (i) collection time; (ii) collection country; (iii) clinical phenotype; (iv) emm type.
FigureS3.pdf
Figure S3: Basic characteristics of mutations detected based on the core genome for the training cohort of 790 GAS strains. (a) Number of detected mutations of the major categories. (b) Distribution of gene-wise mutation counts (left panel) and mutation ratio (right panel). The mean count or mean ratio are shown in the dotted lines. (c) Violin plot presentation of the MAF distribution for each phenotype. (d-i) scatter plots of MAF for pair-wise phenotypes.
FigureS4.pdf
Figure S4: Manhattan plot presentation of the associated variants across the reference genome for all phenotype pairs. The genomic locations are shown for the variants detected using SNP-based GWAS (a), InDel-based GWAS (b), and CDS-based GWAS (c). The genes or gene regions containing the variants are indicated in color.
FigureS5.pdf
Figure S5: Frequency distribution of genomic distance between the non-coding fragments and 5’-ends/3’-end of neighboring genes. The dotted lines show the average distances.
FigureS6.pdf
Figure S6: A pruning and thresholding strategy was applied to SNPs based on LD and functional annotation. (a-j) Frequency distribution of LD block lengths for SNPs with pair-wise correlation r² ≥ 0.8 within a genomic distance of 10 kb. The dotted lines show the average block lengths. (k) Frequency distribution of pair-wise SNP correlation over all phenotype pairs before pruning and thresholding. The dotted line shows the average r².
FigureS7.pdf
Figure S7: Heatmap presentation of the polygenic pattern of GAS phenotypes. The association significance in –log10(p-value) scale is shown between the genes and phenotypes. If one gene contains more than one associated variant, the –log10(p-value) of the most significant variant was used. For compact presentation, a total of 176 genes were selected to present in the figure from the 783 associated genes based on the selection criteria: (i) all genes significant in CDS-based GWAS were selected; (ii) the genes have association relationships in at least five phenotype-pair groups.
FigureS8.pdf
Figure S8: Network presentation of the association between the eight phenotypes and functional categories in the four genomic partitions. (a) intactCDS, (b) trunCDS, (c) 5’-nonCDS, (d) 3’-nonCDS. The weights of the network edges are proportional to –log10P and the node sizes are proportional to the number of genes in each functional category. The positive and negative association relationships were separately considered and indicated in edge color.
FigureS9.pdf
Figure S9: The number of genetic variants used for calculation of explained phenotypic variance and constructing the GPS models for each phenotype pair. The numbers are cumulative by adding additional variant types.

Download PDF

Editor assigned by journal
13 Sep, 2024
Submission checks completed at journal
09 Sep, 2024
First submitted to journal
07 Sep, 2024

You are reading this latest preprint version

compGWAS: a new GWAS tool allows revelation of the genetic architecture and risk stratification for the versatile pathogen Streptococcus pyogenes

Status:

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Methods and Materials

Genomic data preparation and preprocessing

Data cleaning, core genome generation, and mutation detection

Genotype transformation for the new GWAS methods CDS2GWAS and IndelWinScan

GWAS analysis

Variant screening for genomic partitioning and functional analysis

Estimation of variance explained by the genetic variants

Construction of GPS prediction models and performance evaluation

Results

Dataset overview

Mutation detection

The variants in different genomic partitions are involved in distinct biological functions

Genetic architecture of GAS phenotypes is a hybrid model

Pleiotropicity of the associated genes

The genetic variants from three GWAS modalities explained a large proportion of GAS phenotypic variance

Development of a GWAS polygenic score model for phenotypic differentiation of GAS

Discussion

Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1