Disease-drug pairs revealed by computational genomic connectivity mapping on GBA1 deficient, Gaucher disease mice.

We have reported that, in addition to recapitulating the classical human Gaucher disease (GD1) phenotype, deletion of the glucocerebrosidase (GBA1) gene in mice results in the dysfunction of a diverse population of immune cells. Most of immune-related, non-classical features of GD1, including gammopathies and autoimmune diathesis, are resistant to macrophage-directed therapies. This has prompted a search for newer agents for human GD1. Here, we used high-density microarray on splenic and liver cells from affected GBA1(-/-) mice to establish a gene "signature", which was then utilized to interrogate the Broad Institute database, CMAP. Computational connectivity mapping of disease and drug pairs through CMAP revealed several highly enriched, non-null, mimic and anti-mimic hits. Most notably, two compounds with anti-helminthic properties, namely albendazole and oxamniquine, were identified; these are particularly relevant for future testing as the expression of chitinases is enhanced in GD1.

Introduction

Gaucher disease (GD1) is the archetype lysosomal storage disorder resulting from the defective activity of acid β-glucosidase arising from biallelic mutations in the GBA1 gene [1]. Its phenotypic presentation is diverse displaying imperfect correlations with mutations, mainly due to modifier genes [2], [3], [4]. Specifically, certain disease pathologies, such as malignancies [5], gammopathies [6], autoimmune diathesis [7], Parkinson disease [8], and osteoporosis [9] cannot be explained by the lysosomal accumulation of GBA1 lipid substrates in macrophages. Indeed, we have shown that mice lacking the GBA1 gene not only recapitulate the visceral manifestations of GD1, including hepatosplenomegaly, marrow infiltration, and extramedullary hematopoiesis, but also present with widespread dysfunction of unexpected immune cell populations [10]. Despite this new information, we know from clinical studies that most non-classical clinical features of GD are resistant to macrophage-directed enzyme replacement therapy, prompting an exploration into other potential molecules for treatment of GD.

Attempts at drug discovery have traditionally utilized high-throughput screening of small molecule libraries. In contrast, the Connectivity Map (CMAP; http://broad.mit.edu/cmap) is a rapid in silico technique for studying connections between mechanisms of drug action, consequences of genetic perturbations, and molecular aberrations in disease states [11]. Superseding gene profiling, which was used to elucidate biological pathways and reveal cryptic disease subtypes, this new approach employs a disease (or drug) gene signature to query a database comprising gene expression profiles [11]. Non-parametric, rank-based pattern-matching using the Kolmogorov–Smirnov statistic reveals compounds that share mechanisms that overlap those of a given disease or genetic perturbation. These compounds can subsequently be tested for biological activity, potentially as lead molecules.

In the current study, we first developed a query signature from high density microarray profiling of spleen and liver cells obtained from GD1 and wild type mice [10]. We then interrogated CMAP with this pathway-based query signature to yield a hierarchical list of novel compounds that share common mechanisms with GBA1. These compounds will further be tested for their biological activity in vitro and in vivo in GBA1-deficient GD1 mice.

Materials and methods

We used a previously published dataset from the GD1 mouse, in which the GBA1 gene was deleted in cells of the hematopoietic and mesenchymal lineage [10]. GD1 mice exhibited differences in phenotypic severity, and were classified as severe or mild/moderate based on the fold-increase in spleen size compared with control unaffected mice [10]. Fig. 1 shows our overall approach for selecting affected genes, filtering, and connectivity mapping.

Fig. 1

Schematic diagram showing our approach to drug mimic/anti-mimic identification utilizing a combination of techniques, namely, microarray gene profiling, pathway analysis and connectivity mapping.

Gene expression patterns in the liver and spleen were profiled using the GeneChip Mouse Gene 1.0 ST Array (Affymetrix, Santa Clara, CA) and analyzed with the Bioconductor GeneChip Robust Multi-array Average (GC-RMA) package. The dataset is available at the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) (accession number GSE23086).

Genes that satisfied the following criteria were selected as query for Ingenuity Pathway Analysis (IPA) (Ingenuity Systems, Redwood City, CA): (1) their expression pattern followed a progressive change, i.e. severe > mild/moderate > control; and (2) the up- or down-regulated genes in the severe GD1 mouse showed more than a 2-fold change compared with control. After IPA pathway analysis, genes that constituted the most significant network, i.e. having the lowest p-value, were selected. However, not all selected genes were present in the original query as some were merely connecting molecules within the network. We therefore used more stringent filtering criteria for the IPA-analyzed genes, namely that (1) their expression pattern must follow a progressive change, and (2) the up- or down-regulated genes in the mild/moderate GD1 mouse must show a >2-fold change compared with control mice. Note that these criteria are more stringent in that the latter comparisons are being made between control and mildly/moderately affected mice (spleen size <6 multiples of normal), whereas in the initial screen control mice were being compared against severely affected animals (spleen size >6 multiples of normal).

These sets of genes, namely liver up-regulated, liver down-regulated, spleen up-regulated, and spleen down-regulated, were translated to their respective Affymetrix GeneChip Human U133 probe set names and used as the “query signature” for analysis with CMAP (Connectivity Map; build 02; http://broadinstitute.org/cmap/) [11]. Results of the CMAP analysis were ranked first by enrichment and then by percent non-null after discarding hits with p-value > 0.05 and hits with number of instances (n) smaller than 2.

Results and discussion

There are >28,000 probe sets in the Mouse Gene 1.0 ST Array. Using a 2-fold threshold, we found 974 up-regulated and 187 down-regulated genes in the liver, and 546 up-regulated and 499 down-regulated genes in the spleen in the severe GD1 mice compared to the control (Supplementary Tables 1–4). Interrogation of IPA with these gene sets yielded several most significant networks (Supplementary Table 5). The liver up-regulated genes comprised inflammatory response and immune genes, whereas lipid metabolism and fatty acid metabolism genes were down-regulated. In contrast, gene pathways up-regulated in the spleen included infectious disease and severe acute respiratory syndrome genes, whereas those down-regulated comprised genes related to hematological system development, and the function and quantity of lymphocytes. The gene expression pattern that we deciphered is broadly similar to that reported in another mouse model of GD1 [12].

The data are biologically plausible, in view of our recent observation of a plethora of immune cell defects in GD1 mice with the GBA1 gene deleted in hematopoietic stem cells (HSCs) [13]. Notably, we found that the thymus exhibits the earliest and most striking alterations suggesting impaired T-cell maturation, aberrant B-cell recruitment, and enhanced antigen presentation [10], [13]. In contrast, the spleen showed less impressive changes in immune cell composition; these changes were restricted mainly to severely sick mice. Thus, we argue that genetic aberrations in immune gene expression revealed on high density microarray might precede, and perhaps even underlie, the cellular composition changes we note in the spleen. This is a prelude to future studies on human GD1 to evaluate whether query signatures could perhaps even be assigned to predict immune-related manifestations.

As noted above, macrophage-targeted enzyme replacement therapy does not ameliorate, in particular, the immunological and some of the skeletal manifestations of GD1 [14]. We therefore elected to develop spleen- and liver-based query signatures to interrogate CMAP, so as to yield compounds amenable for further testing. For this, we further filtered gene sets contained within the up-regulated down-regulated pathways noted in Supplementary Table 5, specifically to remove genes that were not regulated in the original dataset (Fig. 2, Fig. 3 , Supplementary Table 6). Table 1 shows the number of genes before and after filtering; these gene sets were used as “signatures” to query CMAP.

Fig. 2

Pathway-enriched genes that show up (A) and down (B) regulation in the liver of severe GD1 mice after pathway analysis. Genes were further selected such that they showed at least 2-fold up (C) or down (D) regulation in the liver of mild GD1 mice.

Fig. 3

Pathway-enriched genes that show up (A) and down (B) regulation in the spleen of severe GD1 mice after pathway analysis. Genes were further selected such that they showed at least 2-fold up (C) or down (D) regulation in the spleen of moderate GD1 mice.

Table 1

The number of genes that constituted the most significant network after IPA pathway analysis, and after filtering for ⩾2-fold regulation in the mild/moderate GD1 mouse.

	Pathway analysis		After filtering
	Up	Down	Up	Down
Liver	249	25	22	1
Spleen	22	52	14	35

Table 2 shows the most significant hits after CMAP analysis (see Supplementary Tables 7 and 8 for complete results). The chemical structures of these hits are shown in Supplementary Fig. 1. Of note is that the data is ranked first by enrichment and then by percent non-null hits. This criteria is different from that used previously, wherein connectivity scores were the sole determinants of the CMAP ranking [19]. In line with inflammatory and immune gene regulation, CMAP yielded a host of chemicals utilized in acute and chronic infections. One such high scoring hit, albendazole (enrichment 0.823, p 0.0112, n  = 3), is used as an anti-helminthic agent in cysticercosis. Interestingly, the highest ranking anti-mimic was oxamniquine, again an anti-helminthic agent against schistosomiasis.

Table 2

The 20 most significant hits, mimics and anti-mimics, in GD1 mouse liver and spleen after CMAP analysis. Shown is the mean connectivity score, the number of instances (n), the enrichment score, p-values, specificity score, and percent non-null. The known biological action of each compound is stated. Please refer to Supplementary Fig. 1 for chemical structures.

	CMAP name	Mean	n	Enrichment	p	Specificity	Percent non-null	Biological action
Liver
Mimics	5248896	0.692	2	0.946	0.00547	0.0063	100	–
	5224221	0.693	2	0.881	0.02875	0.2737	100	–
	Albendazole	0.625	3	0.823	0.01112	0.0116	100	Anthelmintic
	Benzethonium chloride	0.599	3	0.779	0.02193	0.0474	100	Antimicrobial
	Alimemazine	0.599	4	0.776	0.00477	0.0121	100	Antipruritic
	Sulconazole	0.455	4	0.719	0.01247	0.0500	75	Antifungal
	Protoveratrine A	0.242	4	0.719	0.01273	0.0074	50	Vasodilator
	8-azaguanine	0.430	4	0.715	0.01331	0.1183	75	Antineoplastic
	Octopamine	0.400	4	0.688	0.01969	0.0000	75	Sympathomimetic
	Saquinavir	0.450	4	0.688	0.01975	0.0074	75	Antiretroviral
Anti-mimics	Oxamniquine	−0.560	4	−0.876	0.00052	0.0000	100	Anthelmintic
	Timolol	−0.542	4	−0.870	0.00058	0.0000	100	Beta-adrenergic receptor blocker
	Yohimbic acid	−0.485	3	−0.842	0.00793	0.0118	100	Aphrodisiac
	Fendiline	−0.432	3	−0.832	0.00953	0.0191	100	Non-selective calcium channel blocker
	Iohexol	−0.536	4	−0.831	0.00147	0.0000	100	Contrast agent
	Ribavirin	−0.433	4	−0.825	0.00177	0.0132	75	Antiviral
	Ketotifen	−0.427	4	−0.824	0.00181	0.0082	75	Antihistamine
	Isoflupredone	−0.427	3	−0.822	0.01126	0.1667	66	Anti-inflammatory corticosteroid
	Methoxamine	−0.422	4	−0.821	0.00195	0.0065	75	A1-adrenergic receptor agonist
	Cefaclor	−0.237	4	−0.817	0.00209	0.0000	50	Cephalosporin antibiotic
Spleen
Mimics	W-13	0.603	2	0.863	0.03811	0.0307	100	Calmodulin antagonist
	Gabexate	0.631	4	0.855	0.00058	0.0000	100	Serine protease inhibitor
	MK-886	0.577	2	0.850	0.04529	0.0435	100	Inhibitor of leukotriene biosynthesis
	Oxolamine	0.410	4	0.793	0.00362	0.0275	75	Cough suppressant
	Triflupromazine	0.482	4	0.762	0.00609	0.0357	75	Antipsychotic
	Acenocoumarol	0.611	5	0.752	0.00228	0.0000	80	Anticoagulant
	Rifampicin	0.503	4	0.728	0.01092	0.0233	75	Bactericidal antibiotic
	Arachidonic acid	0.293	3	0.717	0.04473	0.0513	66	Polyunsaturated omega-6 fatty acid
	Metergoline	0.312	4	0.698	0.01719	0.0812	50	Psychoactive drug
	Meprylcaine	0.380	4	0.697	0.01763	0.0206	50	Local anesthetic
Anti-mimics	Scopoletin	−0.758	2	−0.924	0.01181	0.0144	100	Coumarin
	Skimmianine	−0.786	4	−0.918	0.00010	0.0000	100	Histamine release blocker
	Canavanine	−0.735	3	−0.880	0.00345	0.0052	100	Non-proteinogenic a-amino acid
	16,16-dimethylprostaglandin	−0.649	3	−0.821	0.01134	0.0138	100	Prostaglandin E receptor agonist
	Dequalinium chloride	−0.712	4	−0.805	0.00280	0.0058	100	Antiseptic
	Alfaxalone	−0.598	3	−0.780	0.02189	0.0266	100	Neurosteroid general anesthetic
	Phenindione	−0.619	4	−0.771	0.00557	0.0108	100	Anticoagulant
	Aminoglutethimide	−0.394	3	−0.747	0.03295	0.0238	66	Anti-steroid
	Quinidine	−0.482	3	−0.730	0.04052	0.0301	66	Class I antiarrhythmic agent
	Dexpanthenol	−0.337	4	−0.713	0.01377	0.0286	50	Humectant/emollient/moisturizer

Notable is that the response of the macrophage to the accumulation of major (GL1) and minor (Lyso GL1) lipid substrates arising from GBA1 deficiency has an interesting parallel with infection by chitin-containing parasitic organisms. For example, in addition to a trend towards alternative macrophage polarization [15] and Th2 responses [13], there is an impressive production of chitinases and chitinase-like proteins in GD1 [16], [17]. In fact, there is a massive up-regulation of chitotriosidase in lipid-laden Gaucher cells with a ∼1000-fold elevation in serum levels in human GD1 [16]. Likewise, we have reported a ∼25-fold increase in the expression of chitinase 3-like 1 and chitinase 3-like 3 in the liver and the spleen in GD1 mice [10]. This prompts the future testing of anti-helminthic agents, such as albendazole and oxamniquine, in GD1. A similar example of drug-disease pairing using an inflammatory bowel disease (IBD) gene signature has led to the identification of topiramate, an anti-convulsant, as an effective agent for IBD [18]. Likewise, it has been shown using CMAP that certain HDAC inhibitors have mechanistic overlaps with selective estrogen receptor modulators [19].

In addition to anti-helminths, another notable high-ranking anti-mimic in the liver CMAP output was fendiline (enrichment −0.832, p 0.00953), a calcium channel blocker. The latter class of drugs is already being explored for therapeutic potential in GD1 [20]. Finally, there was no particular pattern to the CMAP mimics or anti-mimics for the spleen; these ranged from a histamine release blocker (skimmianine) to a calcium-calmodulin inhibitor (W-13).

Thus, computational genomic connectivity mapping provides a unique and purposeful strategy for unraveling new actions of drugs currently used for other medical conditions, as well as, new therapeutics for old diseases. Given that the therapeutic armamentarium for GD1 is small, and enzyme replacement therapy does not reverse all GD1 manifestations, it would be worthwhile investigating high-ranking CMAP outputs for biological activity in GBA1 deficient mice.