Complex polymorphisms in endocytosis genes suggest alpha-cyclodextrin as a treatment for breast cancer.

Most breast cancer deaths are caused by metastasis and treatment options beyond radiation and cytotoxic drugs, which have severe side effects, and hormonal treatments, which are or become ineffective for many patients, are urgently needed. This study reanalyzed existing data from three genome-wide association studies (GWAS) using a novel computational biostatistics approach (muGWAS), which had been validated in studies of 600-2000 subjects in epilepsy and autism. MuGWAS jointly analyzes several neighboring single nucleotide polymorphisms while incorporating knowledge about genetics of heritable diseases into the statistical method and about GWAS into the rules for determining adaptive genome-wide significance. Results from three independent GWAS of 1000-2000 subjects each, which were made available under the National Institute of Health's "Up For A Challenge" (U4C) project, not only confirmed cell-cycle control and receptor/AKT signaling, but, for the first time in breast cancer GWAS, also consistently identified many genes involved in endo-/exocytosis (EEC), most of which had already been observed in functional and expression studies of breast cancer. In particular, the findings include genes that translocate (ATP8A1, ATP8B1, ANO4, ABCA1) and metabolize (AGPAT3, AGPAT4, DGKQ, LPPR1) phospholipids entering the phosphatidylinositol cycle, which controls EEC. These novel findings suggest scavenging phospholipids as a novel intervention to control local spread of cancer, packaging of exosomes (which prepare distant microenvironment for organ-specific metastases), and endocytosis of β1 integrins (which are required for spread of metastatic phenotype and mesenchymal migration of tumor cells). Beta-cyclodextrins (βCD) have already been shown to be effective in in vitro and animal studies of breast cancer, but exhibits cholesterol-related ototoxicity. The smaller alpha-cyclodextrins (αCD) also scavenges phospholipids, but cannot fit cholesterol. An in-vitro study presented here confirms hydroxypropyl (HP)-αCD to be twice as effective as HPβCD against migration of human cells of both receptor negative and estrogen-receptor positive breast cancer. If the previous successful animal studies with βCDs are replicated with the safer and more effective αCDs, clinical trials of adjuvant treatment with αCDs are warranted. Ultimately, all breast cancer are expected to benefit from treatment with HPαCD, but women with triple-negative breast cancer (TNBC) will benefit most, because they have fewer treatment options and their cancer advances more aggressively.

Introduction

Breast cancer is the most common cancer in women worldwide.[1] In 2016, 246,660 new U.S. cases were estimated.[2] The highly penetrant, but rare mutations in BRCA1 and BRCA2 point to DNA repair deficiencies as an etiological factor, but explain only 5 to 10 percent of cases. Patients with breast cancer positive for estrogen receptor (ER) or human epidermal growth factor (GF) receptor type 2 (HER2) initially respond well to anti-estrogen or anti-HER2 therapy, respectively, but inevitably become refractory.[3]

As of May, 2016, the deadline for participation in the National Cancer Institutes’ “Up For A Challenge” (U4C) breast cancer challenge, 127 single nucleotide polymorphisms (SNPs) had been associated with breast cancer in women of European ancestry [4] at the conventional fixed s = −log(p) = 7.3 level for genome-wide statistical significance (GWS) [5] (s is used throughout for significance). These SNPs map to 51 genes with known function; all but 16 involved in three known pathways: 27 are associated with nuclear function (DNA repair, transcription, cell-cycle control), six with receptor signaling, ion channels, and mammary gland development (KEGG pathway hsa04915) and two with AKT signaling (hsa04064).[6] The U4C aimed to generate novel testable biological hypotheses (80 FR 32168).

The present evaluation is based on separate analyses of three independent populations of women of European ancestry (see Subjects). Two of the populations (EPIC, PBCS) had never been analyzed individually, because their sample size was deemed insufficient for conventional statistical approaches.

Most breast cancer deaths are not due to the primary tumor, but to metastases, often in the bone, lung, liver, and brain. The genetics results submitted under the U4C implicate dysregulation and dysfunction of endo-/exocytosis (EEC), which is involved in cell migration and invasion, as well as organ targeting, and, thus, suggest overall downregulation of phosphoinositides (PI) as a novel treatment strategy against metastases. The hypothesis that alpha-cyclodextrin (αCD), which scavenges phospholipids, is effective in reducing migration of breast cancer tumor cells was subsequently confirmed in an in vitro study. Taken together, the results suggest (derivatives of) αCD as a potential treatment for carcinomas without the side effects of radiation and cytotoxic drugs or radiation.

Materials and methods

Ethics statement

The study was approved by The Rockefeller University IRB on Aug 24, 2015 (ref# 330390, exempt).

Subjects

This reanalysis is based on data from three GWAS in mostly ER− (including some PR− and/or HER2−) women of European ancestry:

(a) the NHS cases from the Nurses’ Health Study as part of the Cancer Genetic Markers project (CGEM, phs000147/39389-2/GRU, 1145 cases / 1142 controls),[7, 8]

(b) ER− cases from the nested case-control study of the European Prospective Investigation into Cancer (EPIC, phs000812/39395-2/HMB-PU, 511 cases / 500 controls),[9]

(c) ER− cases from the Polish Breast Cancer Case-Control Study (PBCS, phs000812/39397–2, 543 cases / 511 controls),[9]

The EPIC and PBCS studies are part of the Breast and Prostate Cancer Cohort Consortium GWAS (BPC3), which was supported by the National Cancer Institute (NCI) under cooperative agreements U01-CA98233, U01-CA98710, U01-CA98216, and U01-CA98758 and the Intramural Research Program of the NCI, Division of Cancer Epidemiology and Genetics (see https://www.synapse.org/#!Synapse:syn3157598/wiki/232630 for further details).

Statistical methods

In this analysis, conventional single-SNP GWAS (ssGWAS) are complemented with a computational biostatistics approach (muGWAS, GWAS using muStat [10]) that incorporates knowledge about genetics into the method (see Sections 4.3.4 and 4.4.2 in [11]) [12] and knowledge about the nature of GWAS into the decision strategy.[13]

Statistical methods tend to have higher power if they are based on more realistic assumptions, which, in biology, tend to be weak. In contrast, methods based on stronger assumptions, such as additivity of allelic effects and independence of SNPs within an linkage disequilibrium (LD) block (LDB), may generate more significant results when errors happen to fulfill these assumptions than for true effects. With millions of test statistics calculated, even a rare false positive result due to model-misspecification (1/10,000 tests, say), may result in the 100 most significant results all being false positives. U-statistics for multivariate data in GWAS (muGWAS) rely only on weak, realistic assumptions, but require large amounts of memory and GPU enabled cloud instances, which became available only after 2001 and 2009, respectively.

After excluding non-informative or low-quality SNPs and SNPs in high LD with an immediate neighbor [14] (20–25%) to avoid loss of power when including irrelevant SNPs [15], an initial traditional ssGWAS was performed, using the u-test for univariate data.[16-18] The same data was then analyzed using a u-test for genetically structured multivariate data.[12] U-statistics avoid model-misspecification biases by replacing linear/logistic [19] with non-parametric kernels.[15]

Below, we describe the assumptions about genetics and GWAS that are implemented in the statistical method and decision strategy and refer to published empirical validation of this approach.

1.1 Heterodominance

A particular SNP is not assumed to be either recessive (aA = aa), additive (aA = (aa+AA)/2), or dominant (aA = AA), but merely monotonic (aa < aA < AA). Accordingly, the information contributed by a particular SNP is represented as a matrix detailing for each of the n×n pairs of n subjects whether the genetic risk carried by the row subject is lower “<”, the same “=“, or higher “>” than the risk of column subject, or unknown (“?”) in case of missing data in one or both of the subjects. Below, the possible genetic risk constellations (left) are compared to models with different degrees of dominance (right). While the left matrix is similar to the matrix for dominant effects (all non-zero elements are ±2), the (logical) inequalities are not (numerically) equivalent. In effect, the single-SNP results based on the adaptive u-scores for aa, aA, and AA are similar to results from the Cochran-Armitage test for additive co-dominance, [20, 21] which uses fixed scores 0, 1, and 2.

1.2 LD-structure

A basic assumption underlying GWAS, in general, is that a disease locus should be in LD with both neighboring SNPs (unless they are separated by a recombination hotspot). Hence, the information from two neighboring SNPs is not numerically ADD-ed, but logically AND-ed using the function Ʌ

As muStat allows variables to be correlated, other SNPs within an LDB may be in LD, too, yet there is no formal representation of more distant LD. Non-informative SNPs added between LDBs prevent intervals from spanning LDBs.

1.3 Cis-epistasis, including compound-heterozygosity

To account for interactions between functional polymorphisms,[22] a natural extension of Ʌ is then used to combine information from corresponding elements of the n×n matrices containing information about neighboring pairs. Assuming, without loss of generality, the case of only four SNPs within in the same LDB, the aggregated diplotype information for one pair of subjects is

which can be one of the following (invariant to permutaitons π)

From the above inequality, the results typically differ when SNPs from the same tag sets appear in different permutations, which increases the resolution over methods assuming commutativity.

1.4 Test statistic

From the resulting n×n matrix W (say), one calculates each subject’s risk score u (−n < u < n) as the number of subjects having lower risk, minus the number of subjects having higher risk, i.e., #(w = “<”)j − # (w = “>”)j. These scores are then compared between cases and controls using a standard linear score test.[23]

1.5 Regularization

Since it is unknown a priori, whether a minor allele is dangerous, irrelevant, or protective, all combinations of (−1, 0, +1) “polarities” are applied to the SNPs S, …, S, resulting in many highly dependent test statistics being calculated for the diplotypes surrounding a given SNP. The test statistic chosen is the one that has the highest u(−log(p), IC ) score, where the information content (IC) is the proportion of pairwise orderings in W that can be decided (≠“?”) for a given choice of polarities. This approach avoids over-fitting (highly significant results based on a small subset of unusual subjects) without the need to choose arbitrary regularization cut-offs.[24]

2.1. Adaptive genome-wide significance

The traditional p-value cut-off of s = 7.3 for GWS has been widely criticized as overly conservative,[25, 26] yet few alternatives have been formally derived. Here, we replace a fixed cut-off for GWS with an empirical [22] adaptive (study-specific) cut-off (aGWS) that automatically accounts for the specifics of the population studied, the chip used, differences in minor allele frequency (MAF,) and GWAS being non-randomized.[13] As previously discussed,[13] the expected distribution in a ssGWAS QR plot is a mixture of univariate distributions whose carriers vary by MAF, because the most significant result possible depends on MAF when outcomes are bounded (allele counts 0, 1, 2). Hence, it is a convex curve, rather than a straight line;[13] see, for instance, CGEM chromosomes 14–17, 19, and 22 (Fig B in S1 File). In a whole genome (WG) plot, this curvature may not be apparent (see Fig 1, below), when some chromosomes’ QR curves are concave because of true association, which is expected in a familial disease or with systematic unrelated differences between non-randomized populations. Hence, an apparently straight line in a WG plot may be due to concave curves in chromosomes with true positives and convex curves in others canceling each other out. With muGWAS, where many dependent tests are performed at overlapping window positions, the expected QR curve (see Fig C in S1 File) may be even more convex. The expected distribution curve is estimated from the 50% of chromosomes with the fewest outliers rising above a convex fit.[13] The empirical adaptive (study-specific) aGWS cut-off is the median apex (highest point) of a convex curve fitted against these chromosomes’ QR plot.

Fig 1

GWAS quantile-rank (QR) plots.

Left: ssGWAS, right: muGWAS (each point represents the most significant result among all diplotypes centered at the same SNP) Results are ranked by significance (bottom). For the most significant results and other results of interest, the location of SNPs to genes is shown in. Upper curve (red): convex fit against points; dashed extension: projection; lower curve (blue): population-specific expectation. Vertical lines between curves connect the highest s-values (−log10 p) of a gene (dot) with its expected value for genes with known function. Light gray vertical lines indicate genes omitted from the list because of low reliability (either low µIC or reliance on a single SNP), respectively. Genes to the right of the vertical dark line are above the aGWS cut-off. See Fig A in S1 File for Manhattan plots. The horizontal solid line at highest point at the end of the expected curve indicates the estimate for adjusted GWS (aGWS). All points above the horizontal line (and genes to the right of the vertical blue line) are “significant” at the aGWS level.

2.2. Replication

Complex diseases may involve different SNPs in high LD with causal variants across populations,[27] epistasis between several SNPs per locus, several loci per gene, and several genes per function, with risk factors differing across populations (see above). Hence, we will consider SNPs within a locus, loci within a gene, and genes with a direct mechanistic relationship (paralogs, binding partners,…) for replication. [22, 28] Results are considered “replicated” if supportive results are significant at the aGWS/2 level.

Validation

The above approaches have been validated in two published analyses, where previous analyses using ssGWAS and fixed GWS also had identified not more than a few apparently unrelated SNPs.

In epilepsy,[12] muGWAS confirmed the Ras pathway and known drug targets (ion channels, IL1B). In that analysis, muGWAS was also compared with a parametric analogue, logistic regression with interaction terms for neighboring SNPs (lrGWAS). muGWAS produced fewer apparent false positives (isolated highly significant results far away from coding regions) (Suppl. Fig 2 in [12]) and higher sensitivity for genes downstream of Ras, which are involved in more complex cis-epistatic interactions, (Fig 3, blue, in [12]) than ion channels, which were also implicated by lrGWAS (see Fig 3, red, in [12]).

In autism,[13] muGWAS identified sets of mechanistically related genetic risk factors for mutism in autism (independently confirmed in functional studies [29] and a pathway network analysis [30]). In, [13], adaptive GWS was validated against three analyses with randomly permutated phenotypes. Only one gene (DMD, not aGWS) appeared in one of the other analyses (also not aGWS). Moreover, there is no noticeable overlap between aGWS genes between breast cancer and either mutism [13] or epilepsy (Suppl. Fig 7 in [13]), while there is considerable functional overlap between mutism in autism and epilepsy, as expected.

In vitro Assay

A 24-well plate (CBA-120, Cell BioLabs Inc., San Diego, CA) with CytoSelect Wound Healing Inserts was warmed up at room temperature for 10min. A cell suspension used contained 0.5–1.0×106 cells/ml in media containing 10% fetal bovine serum (FBS) was prepared and 1mL of this suspension was added to each well. Cells were then incubated for 12h, after which time the insert was removed and cells were washed with new media to remove dead cells and debris. FBS with/without CDs (Sigma-Aldridge, St. Louis, MO) was added to start the wound healing process. Cells were incubated for 2h, washed with PBS, fresh control media was added, and cells were incubated for another 12h. After removing the fixation solution, 400µL of Cell Stain Solution were added to each well and incubated for 15min at room temperature, after which stained wells were washed thrice with deionized water and left to dry at room temperature. Cells that migrated into the wounded area or protruded from the border of the wound were visualized and photographed under an inverted microscope to determine migrated cell surface area.

https://www.cellbiolabs.com/sites/default/files/CBA-120-wound-healing-assay.pdf

Results

Additional ssGWAS CGEM results complement known breast cancer risk factors

The original CGEM analysis had identified two SNPs (rs1219648: s = 5.49, rs2420946: 5.46) in the fibroblast GF receptor FGFR2 Entrez Gene 2263,[7] which affects mammary epithelial cell growth and migration,[31] and a SNP (rs10510126: 6.25, >1 MB apart from FGFR2) which was subsequently located to a long variant of the mitotic checkpoint protein BUB3 9184. These two genes are also the only genes in the present analysis with SNPs above the diagonal in the summary ssGWAS quantile-rank (QR, often: QQ) plot (Fig 1 left), although the QR plots of several individual chromosomes show association in chromosomes 4 (the SNCA–MMRN1 22915 region), 5 (breast cancer associated transcript BRCAT54 100506674, non-coding), 6 (PARK2 5071, the Parkinson’s disease [PD] ubiquitin ligase Parkin), and 9 (LPPR1 54886, phospholipid phosphatase-related 1) (Fig B in S1 File).

In the present analysis, a total of 22 genes and BRCAT49 [32] reached aGWS in CGEM (Fig 1, left, blue). A total of 21, 11, and 24 genes with known function or relation to breast cancer exceeded muGWAS aGWS in CGEM, EPIC, and PBCS, respectively.

Novel ssGWAS aGWS results in EPIC and PBCS complement CGEM results

In EPIC, the two most significant SNPs (rs4791889: 5.66 and rs9596958: 5.42) are located 4.5 kB upstream of the chromodomain helicase DNA binding protein CHD3 1107 and the transcription factor (TF) SOHLH2 54937, respectively (Table B in S1 File and Fig E in S1 File).

In PBCS, the two most significant SNPs (rs2297075: 5.83, rs943628: 5.55, 100 kB apart) are located in DOCK8 81704, a guanine nucleotide exchange factor for Rac1, which drives mesenchymal cell movement.[33] Significance of FGFR2 relies on the two previously reported and a third SNP (rs11200014) within intron 2.[34] Significance in BUB3 is driven by three SNPs in high LD (rs10510126, rs17663978, rs7916600, spanning 30 kB). These findings are consistent with significance of the top five SNPs in ssGWAS depending on a single polymorphism each. Lack of evidence in EPIC and PBCS (Table B in S1 File) is consistent with different variations developing in divergent European populations.

muGWAS aGWS results are cross-validated across CGEM, EPIC, and PBCS

In CGEM, the top gene was the phospholipid/diacylglycerol (DAG)-dependent protein kinase PRKCQ 5588 (chr10: 6,540,724–6,573,883), which induces cell migration and invasion.[35, 36] The same SNP (rs661891) was also implicated in EPIC. The three most significant SNPs and the most significant regions in muGWAS were all located within the same 34kB LDB. The second most significant gene was a long EST of the transient receptor potential cation channel TRPM3 80036, which controls oncogenic autophagy in renal cell carcinoma,[37] supported by a part of the promoter region of the shorter main form in PBCS. BUB3 was also significant in muGWAS, followed by the endo-/lysosomal receptor SCARB2 950 and the nuclear RNA polymerase subunit POLR1A 25885 (rs10779967).

In EPIC, the top gene in muGWAS (as in ssGWAS), was the TF SOHLH2, followed by AGPAT3 56894 (rs8132053 in CGEM and EPIC), whose paralog AGPAT4 56895 is included in Fig 1 (4.94, right panel, CGEM). CELF2 10659, an RNA binding protein, and STARD13 10948, a breast cancer tumor suppressor that regulates cell migration and invasion [38] also reached aGWS. CHD 364663 depends entirely on SNP rs4791889 (see Statistical Methods, 2.2. Replication, for replication criteria).

In PBCS, the top gene in both muGWAS and ssGWAS was DOCK8 81704, followed by the nuclear receptor corepressor NCOR2 9612, which has been implicated in tamoxifen resistance in breast cancer.[39, 40] CACNA1C 775 (3rd) is highly up-regulated in breast cancer.[41] The multiple epidermal GF-like domains protein 11 (MEGF11 84465, 4th), like MEGF10 84466 an ortholog of C. elegans Ced-1 and the Drosophila draper, had been implicated in colorectal cancer. [42]

Both CGEM and EPIC have a significant P-type ATPase, which import phosphatidylserine (PS, ATP8B1 5205) and phosphatidylcholine (PC, ATP8A1 10396), respectively, the substrates for phospholipase D (PLD) to produce phosphatidic acid (PA) for the synthesis of phosphatidylinositol (PI).[43] BMP7 655 (ss: 4.24) and its receptor BMPR1B 658 (ss: 4.47) are significant in EPIC and CGEM, respectively, and BMP signaling is known to regulates mitotic checkpoint protein levels in human breast cancer cells, including levels of BUB3 (see above).[44] DGKQ 1609 (rs2290405) which converts DAG into PA, was replicated in CGEM and PBCS, while LPPR1 54886, which is involved in the conversion of PA into PI was replicated in CGEM and EPIC.

As expected in samples from the general population, the known risk factors for rare early-onset breast cancer (BRCA1/2 672/675, HER2 2064, RB1 5925) do not show association and many receptor-related genes are absent in ER− populations. Except for the genes with highest significance in ssGWAS (BUB3 in CGEM, SOHL2 in EPIC, and DOCK8 in PBCS), all of the aGWS genes in muGWAS have support in least one of the other two populations (2nd block of Table B in S1 File). This observation is consistent with muGWAS identifying primarily old cis-epistatic variations, rather than de novo mutations favored by ssGWAS. Table B of S1 File gives an overview about the significance and replication of the genes identified and supportive evidence in the literature.

muGWAS results confirm known disease pathways in breast cancer

Consistent with the published results in the NHGRI-EBI catalog, a total of 16, 15, and 18 genes above aGWS in CGEM, EPIC, and PBCS, respectively, are involved in the three known disease pathways, such as membrane-associated receptor signaling (G protein–coupled receptors [GPCR], Fc receptors [FcR], hemagglutinin [HA], receptor tyrosine kinases [RTK], or ion channels), MAP kinases, and in nuclear proteins involved in cell cycle control, transcription, or splicing in breast cancer (Table 1).

Table 1

Breast cancer genes associated with pathways by study.

Within each study (major columns), genes are grouped by function. Mbrn: membrane-associated (GPCR, FcR, HA, RTK, Ion channels), Ncls: nuclear (cell cycle control, transcription, splicing), MPK: MAP kinases, PI/EC: PI cycle/EEC,.Othr: other. Within each block, muGWAS genes (Fig 1) are sorted from top by s-value (s6). s-values above aGWS (CGEM: 5.29, EPIC: 5.71, PBCS: 5.13) are shown in bold. Genes above aGWS in ssGWAS only (CGEM: 4.03, EPIC: 4.00, PBCS: 3.84) are sorted from bottom up (s1); ssGWAS results for genes also implicated in muGWAS are shown next to the muGWAS results. See Table A of S1 File for Entrez Gene identifiers and Table B of S1 File for replication across populations, which is indicated in bold names.

	s6	s1	CGEM		s6	s1	EPIC		s6	s1	PBCS
Mbrn	6.26		TRPM3	Mbrn	6.37		PCSK5	Mbrn	7.74	5.83	DOCK8
	6.00	5.57	FGFR2		5.57		TAS2R1		6.59	4.95	CACNA1C
	5.89		GNG7		5.23		NEFL/M		6.22	4.82	GPC6
	5.85		SYK		5.21		FLT3		5.97	5.80	RGS3
	5.79		CALN1		5.16		EPHB1		5.63		SSTR4
	5.44	4.63	HAS2		5.14	4.08	GRIA1		5.43		EDN1
	5.38		KCND3		5.13		RASD2		5.35	4.86	TRPC3
	5.14		SHANK2		4.93		ACAN		5.22	4.13	HAPLN3
	5.12		LYZL1		4.90		FKBP1A		5.21	5.08	PTPRG
	4.98		LPHN3		4.83	4.24	BMP7		4.95		SDK1
	4.82	4.06	NTSR1		4.64	4.26	GPHN		4.38	4.38	SYNDIG1
	4.49	4.07	MRGPRE		3.68	4.20	TNFRSF10A
					3.57	4.16	PNOC	Ncls	6.63		NCOR2
Ncls	6.20	6.20	BUB3						5.39		ANAPC13
	5.95	4.89	POLR1A	Ncls	8.58	5.42	SOHLH2		5.22	4.46	RBM23
	5.81		MSL3		6.51	4.11	CELF2		5.13	4.39	CDK18
	5.49		TCF15		6.48	4.56	STARD13		5.05		PAX2
	5.42		NUP54		6.13	5.66	CHD3		5.02		GLIS2
	5.22		MACROD2		5.82	4.34	CDKAL1		4.98		ZFAT
	5.08		CTNNBL1		5.49	4.05	SIX2		4.92		HDAC4
	5.06		SAMHD1		5.24		E2F3		4.91		MDFI
	5.03		DMRTB1		5.21		NCCRP1		4.91		TFAP2A
	5.02		ZMAT4		5.08		RARB		5.00	3.83	ISM1
	4.98		SCMH1		5.05		HMGXB4		4.12	4.12	A2BP1
	4.97		HNF1B		5.02		LZTS1		4.61	4.48	FHIT
	4.97		SMARCAL1		3.90	4.18	SFRS8
	4.94		ZBTB20					MPK			—
	4.29	4.29	FBXO38	MPK	5.85	4.36	TRAPPC9
					5.26		PRKCQ	PI/EC	6.39	4.07	MEGF11
MPK	6.70		PRKCQ		5.03	4.32	PRKAG2		6.14		PTENP1
									5.32	4.27	EEA1
PI/EC	6.02		SCARB2	PI/EC	6.59	4.73	AGPAT3		5.17	4.42	ANXA4
	5.58		ATP8B1		6.26	4.87	ATP8A1		4.91	4.41	UNC13C
	5.48	4.65	LPPR1		5.88	4.19	TNS1		4.84	4.14	RAB32
	5.36		SYT17		5.40		GLB1		4.59	4.07	N4BP3
	5.30		SYNJ2		5.33		ASTN2
	5.06		MEGF11		5.27		CHMP7	Othr	6.06	4.60	CDH4
	5.04		NLRP4		5.05		VAV3		6.05		CLLU1
	4.99		ABCA1		4.90		SDCBP2		5.70		TMCC3
	4.94		AGPAT4		4.22	4.22	RAPGEF4		5.34	4.15	COL11A1
	4.59	4.11	STXBP1						5.33		PRH1
	4.97	4.18	SCMH1	Othr	6.33	4.68	mirLET7B		5.21		RNASE11
	4.21	4.21	ANO4		5.50		TMEM132C/D		5.14		CADM2
	4.86	4.24	HNF1B		5.08		COL3A1		5.12		SLC13A3
					5.06		SLC43A3		5.11		KIF25
Othr	5.86		VWA3B		5.04		VWC2L		5.07		MMD2
	5.38	4.46	BRCAT49		5.00		LARGE		5.05		TRIM36
	5.37		MDGA1		5.00		LSMD1		5.02		CYB5RL
	5.00		ALPP		4.98		IGSF9B		4.91		HLA-C
	4.04	4.04	MRPL35		4.88		KPNA3		4.98	4.19	ZNF98
	4.28	4.19	DAP		4.81	4.06	USP44		4.74	4.28	SLC25A26
	4.24	4.19	ABO		5.15	4.12	CFAP36		4.95	4.29	TRMT11
	5.13	4.24	AK5		5.35	4.13	SPATA19		4.42	4.41	CHIT1
	4.89	4.47	BMPR1B		4.20	4.20	MICAL2		4.99	4.45	BCL10
	4.52	4.47	PTCD3		4.77	4.77	OR52K2
	4.95	4.81	BRCAT54
	5.09	4.95	MMRN1

muGWAS results highlight Endo-/Exocytosis (EEC) as a pathway in breast cancer

The cell’s major fibronectin-binding integrin (α5β1) is key to survival and migration of tumor cells.[45] Results of various expression and functional studies have pointed to EEC of β1 integrins as a functional component of “derailed endocytosis” in cancers, including breast cancer (Fig 2).[46-48].

Fig 2

EEC of β1 Integrin underlying mesenchymal tumor cell migration and invasion.

Cell migration necessitates trafficking of β1 integrin, whose internalization is controlled by dynamin. Both clathrin- and caveolin 1 (CAV1)-coated domains of the plasma membrane are involved. Once in early endosomes (EE), integrins may be sorted for degradation in lysosomes, recycled to the plasma membrane through a RAB4-dependent route, or transported to the recycling endosome (RE). Recycling from the RE requires Rab11 family members, such as RAB25 which is often aberrantly expressed in human tumors, including luminal B breast cancer.[49] (adopted from [46–48]).

Among the 15 GWS genes not associated with known pathways in the NHGRI-EBI catalog (excluding the ambiguous locus between MDM4 4194 and PIK3C2B 5287), only four are involved in EEC (PDE4D, SNX32, STXBP4, DNAJC1, marked with “*” in Table A of S1 File), all from ssGWAS of a combined analysis of nine studies,[50] which included the three studies analyzed separately here. A String() pathway analysis of the subset of aGWS genes that are not part of the above three pathways identified two clusters related to EEC (see Fig 3):

Fig 3

Endo-/exocytosis pathway.

Pink: genes identified in this analysis, most of which have been implicated in breast cancer previously (Table A of S1 File), by stage of EEC: Formation of clathrin-coated vesicles, E3 ubiquitination, separation of inactive integrin (fast recycling) from active integrins (slow recycling), sorting between secretory, lysosomal, and (slow) recycling pathway, and lysosomal degradation. Red and green genes are known breast cancer promoters and suppressors, respectively (Table C of S1 File). Clathrin-mediated endocytosis (CME) begins with co-assembly of the heterotetrameric adaptor complex AP-2 with clathrin at PI(4,5)P2-rich plasma membrane sites. AP-2 in its open conformation recruits clathrin and additional endocytic proteins, many of which also bind to PI(4,5)P2. Maturation of the clathrin-coated pit (CCP) may be accompanied by SHIP-2-mediated dephosphorylation of PI(4,5)P2 to PI(4)P. Synthesis of PI(3,4)P2 is required for assembly of the PX-BAR domain protein SNX9 at constricting CCPs and may occur in parallel with PI(4,5)P2 hydrolysis to PI(4)P via synaptojanin, thereby facilitating auxilin-dependent vesicle uncoating by the clathrin-dependent recruitment and activation of PI3KC2α, a class II PI3-kinase. PI(3,4)P2 may finally be converted to PI(3)P en route to endosomes by the 4-phosphatases INPP4A/B, effectors of the endosomal GTPase Rab5. Adapted from [51]. In the early endosome, β1 integrins are sorted. Inactive integrins undergo fast “short loop” recycling; active integrins go to the late endosome / multivesicular body for slow “long group” recycling (RAB11), lysosomal degeneration (unless rescued by RAB25/CLIC3), or secretion via the trans-Golgi-network mediated by RAB9. Fast recycling of epidermal GF receptor drives proliferation,[52] so one would expect gain-of-function mutations in the upper part of the Figure. Lysosomal and synaptic vesicle exocytosis share many similarities. Endolysosome-localized PIs may regulate lysosomal trafficking in early onset lysosomal storage diseases.[53] and, particularly in ageing, insufficient lysosomal degradation contributes to Alzheimer’s disease (PSEN1, GLB1), Parkinson’s disease (ATP13A2),[54] atherosclerosis (vATPase),[55, 56] and type 2 diabetes (GLB1, HEXA).[57], (derived, in part, from KEGG pathways hsa04144, hsa04721, hsa00531, and hsa04142).

EEC Function: PARK2, PTEN (from PTENP1), SYNJ2, STXBP1, UNC13 (consistent with previous functional studies, see Table C of S1 File)

EEC Regulation: AGPAT3 and DGKQ (Fig D of S1 File).

muGWAS identified genes causing dysfunction of EEC, a known BC risk factor

Further String subset analyses and a literature review by the authors identified additional aGWS genes as related to EEC-related KEGG pathways (genes in parenthesis replaced by a related gene with known function in String). They include endocytosis (hsa04144): DNM1 (from MEGF11), EEA1, PDE4D, SNX32, NEDD4 (from N4BP3) (FDR = .018) and synaptic vesicle cycle (hsa04721): STXBP1, UNC13C, VAMP2; (FDR = .0001).

Fig 4 integrates the genes identified in the present GWAS analysis (pink, see Table A of S1 File for details) with results from expression and functional studies of β1 integrin EEC in breast cancer (see Table C of S1 File for details).

Fig 4

Functional relation of the PI/EC genes.

PI is synthesized from myo-inositol (imported by HMIT) and PA (via CDP-DAG) which can be synthesized from lysophosphatic acid (LPA), PC, or PS, or salvaged from IP3 and DAG. It can also be synthesized from 1-acyl GPI. Arrows: PIs are phosphorylated at a 3-, 4-, or 5- position by PI-kinases (left to right) and hydrolyzed by phosphatases (right-to-left). Genes associated with breast cancer in this GWAS are highlighted in pink (bold: aGWS). See Table 1 for other box colors. Colored arrows in the center indicate the sequence of PIs involved in EEC (Fig 3). Percent values indicate approximate proportion of phospholipids.[58].

muGWAS identifies PI cycle dysregulation as novel breast cancer risk factor

In relation to EEC regulation, both CGEM and EPIC identified a phospholipid-translocating ATPase, ATP8B1 (PE) and ATP8A1 (PS), respectively. AGPAT3 is the second most significant gene in EPIC (mu: 6.59, ss: 4.73); AGPAT4 is among the supportive genes in CGEM (Fig 1, mu: 4.94). Both acyltransferases transform LPA into PA. CGEM also identified the scramblase ANO4 121601 (ss: 4.21), a PS exporter, and the plasma membrane PC/PS efflux pump ABCA1 19 (mu: 4.99). For (ATP8A1, ATP8B1, ANO4, ABCA1), String identified functional enrichment in

GO:0097035 (biol. process) Regulation of membrane lipid distribution: FDR = 0.012

GO:0015914 (biol. process) phospholipid transport: 0.0407

GO:0005548 (mol. function) phospholipid transporter activity: 0.00968

As shown in Fig 4 (upper left corner), 8 (including 6 aGWS) genes are involved in providing the PI cycle with its substrate, PI (and the MAPK signaling pathway with PA.(hsa04072)).

Results for EEC regulation and function are consistent across populations

All three populations show aGWS association with EEC genes (CGEM: 4 in ssGWAS only / 4 in muGWAS only / one in both; EPIC: 1/0/3; PBCS: 3/1/3). Most are validated in at least one of the other two populations, either by the same SNP involved (AGPAT3, DGKQ), the same region (SYNJ2, PDE4D), the same gene (see), or a functionally related gene (AGPAT3/AGPAT4, LPPR1/DGKQ, ATP8A1/ATP8B1, STXBP1/UNC13C, TNS1/PTENP1, see Table B of S1 File for details).

PI supply into the PI cycle as a drug target in breast cancer

After loss-of-function in PTEN and gain-of-function in PI3K suggested a mechanism for upregulation of PI(3,4,5)P3 in cancer, blocking PI3K with Wortmannin [59] or related drugs [60] were considered for treatment of cancers, including breast cancer. Upregulation in PI(3,4)P2 (gain-of-function in SYNJ1/2 or INPPL1 [61]) and PI(3)P (gain-of-function in INPP4B),[62] have also been associated with breast cancer. Recently, components to lower PI(3,4)P2 by inhibiting SYNJ2 have been identified.[63]

Targeting individual phosphotransferases is unlikely to succeed given the robustness of the PI cycle.[59] All PIs regulating EEC, except for the evolutionarily recent MTMR1 link (Fig 4), are regulated by both three kinases and three groups of phosphatases. Given the plethora of genes involved in EEC (Fig 3) identifying the appropriate set of phosphotransferase for a given patient to interfere with endocytosis or to correct for functional deficits in exocytosis may be impractical.

Regulating EEC by controlling the availability of phospholipids, however, while leaving functional interactions within the PI cycle intact, may be feasible. In fact, adding of either exogenous PS or PE led to an enhancement of endocytosis.[64] As EEC is an essential and highly conserved mechanism for tissue morphogenesis [65, 66] and neuronal migration,[67-69] loss-of-function mutations would likely terminate embryonal development. Accordingly, the overall effect of the variations identified (Table C of S1 File) is likely gain-of-function.

HPaCD is more effective than HPbCD against migration of breast cancer cells

In 2014, it was reported that the benefit attributed to the neurosteroid allopregnanolone in the treatment of Niemann-Pick type C (NPC) disease was due to the excipient 2-hydroxypropyl-beta-cyclodextrin (HPβCD). Cyclodextrins are hydrophilic rings of ≥6 starch molecules (Fig 5). The lipophilic cavity can transport lipid drugs, such as allopregnanolone. Empty CDs, at therapeutic doses, form a pool in the aqueous phase into which, in the case of βCDs, cellular cholesterol is extracted,[70] the mechanism of action in NPC.[71]

Fig 5

Structure of cyclodextrins.

Cyclodextrins are toroids formed of six (n = 4, αCD), seven (n = 5, βCD), or eight (n = 6, γCD) starch molecules. The cavity is lipophilic, while the surface is hydrophilic.

Given the focus on cholesterol in NPC, it has often been overlooked that βCDs also scavenge phospholipids. The above GWAS results (Table 1) suggested defects in phospholipid, rather than cholesterol function. Hence, the efficacy of HPβCD in breast cancer might be due to its ability to scavenge phospholipids.

HPβCD is known to inhibit migration of human MDA-MB 231 breast cancer cells.(Fig 3B in [72]) [73, 74] To determine whether inhibition of migration is caused by HPβCD depleting cholesterol, as assumed previously, or by it depleting phospholipids, as implicated by the novel genetics results, the published activity from wound healing experiments comparing HPβCD against control was replicated, and complemented with novel activity results comparing HPαCD against control,,both in MDA-MB 231 (ER–) and MCF-7 (ER+) human breast epithelial cell lines.

From Fig 6, 1mM HPαCD is more effective than 2mM HPβCD against migration of ER− and ER+ tumor cells (p = .0252) while more than 10× less toxic,[76] Hence, the effect previously seen with HPβCD is, in fact, likely the effect of it scavenging phospholipids, rather than cholesterol.

Fig 6

Wound healing by cyclodextrins in breast cancer cell lines.

Cells were grown in triplicates for 12h and incubated with either of the CDs for 2h at the concentration indicated (0–4mM), before a 0.9mm wide gap was opened and cells were allowed to migrate into the “wound” for 12h. HPβCD is more than 10× as toxic as HPαCD, which at <100mM does not affect epithelial cell viability.[75, 76] Dashed horizontal line indicates inhibition of wound healing in HPαCD at 1 and 4 mM respectively. ANOVA results: indep: HPαCD vs HPβCD (fixed) block: MCF-7/MDA-MB-231 (fixed) dep: %change in wound healing 1mM α vs 1mM β, p = .0001 ***1mM α vs 2mM β, p = .0252 *4mM α vs 4mM β, p = .0442 *.

Discussion

Our analysis confirmed previous GWAS, which pointed to receptor/AKT signaling and nuclear functions as critical components in breast cancer etiology. The present results from a reanalysis of data found previously inconclusive provides the first GWAS evidence for the contribution of EEC dysfunction and novel evidence for overstimulation of EEC in mesenchymal tumor cell migration and invasion. These findings, confirmed by an in vitro study on the activity of HPαCD vs HPβCD against breast cancer cell migration, suggest the novel hypothesis that reducing the influx of phospholipids, including PS, PC, and lysophosphatidylcholine (LPC), into the PI cycle via HPαCD could provide an urgently needed treatment option for women with breast cancer.

Replication and complementation of previously identified genes

A previous analysis of the CGEM data reported only two genes, FGFR2 and BUB3, as risk factors for breast cancer. The EPIC and PBCS data have been published only as part of four meta-analysis, which also included CGEM. Among ER− cases, the first meta-analysis [9] confirmed two SNPs each in BABAM1 (7.31) (a nuclear BRCA1 complex component), PTHLH (12.8) (which regulates epithelial-mesenchymal interactions during the formation of mammary glands), and the ER ESR1 (9.6). Our findings of BMP7 (EPIC) and BMPRT1B (CGEM) are consistent with the previous finding of PTHLH, which forms a nuclear complex with BMP4. The second meta-analysis,[77] pointed to the PIK3C2B-MDM4 region (11.68), LGR6 (7.85) (a GPCR), and FTO (7.40) (a regulator of nuclear mRNA splicing). Hence, ssGWAS in all three populations point to receptor/AKT signaling and nuclear processes, although the individual genes differ.

Three of the four EEC genes identified in previous ssGWAS;[50] were confirmed in muGWAS at aGWS/2 (CGEM: 2.56 / EPIC: 2.86 / PBCS: 2.57, Table A of S1 File) in regions in LD (r2):[78] PDE4D (rs1353747, 4.56/4.46/2.84, r2 ≤ .213); SNX32 (rs3903072, 2.92/—/—, r2 ≤ .482, rs7114014); STXBP4 (rs6504950 2.85/—/—, r2 ≤ .238); DNAJC1 (rs11814448,—/—/—).

The EEC genes identified in here (with the exception of AGPAT3/4, ASTN2, and EEA1), have previously been shown to be associated with breast cancer in gene expression and functional studies (Table A of S1 File).

A third meta-analysis [79] based the above three and eleven other U4C data sets,[80] identified five novel breast cancer genes, three with nuclear function (RCCD1, ANKL1, DHODH [81]); ACAP1 and LRRC25 were hypothesized to be involved in cell proliferation (activating Arf6 protein) and inflammatory response (activating hematopoietic cells), respectively, [79] In fact, both genes are can be related to EEC/PI in metastases: ACAP1 (Fig 3, top right) regulates recycling of integrin β1 during cell migration [82]; LRRC25, which regulates development of neutrophils needed for metastases,[83] carries a PI3K interaction motive.[84]

A fourth single-SNP meta-analysis of 68 studies (including the three studies separately analysed here) with a total of 227,000 subjects [85] identified “65 new breast cancer risk loci” to be “incorported into risk prediction models”, with “exocytosis” as the second-most significant “theme” ([85], Suppl Tab. 24).

Computational biostatistics approach to genetic data

The analysis approach,[12] used here integrates genetics concepts into the statistical method, rather than considering them during visual inspection of p-values calculated one SNP at a time and correlations among SNPs within an LDB. In particular, muGWAS avoids assumptions about the degree of dominance, reflects that both SNPs next to a disease locus should be in LD (unless they are separated by a recombination hotspot), increases resolution within LDBs (by distinguishing between members of the same tag sets being in a different order), integrates information from different disease loci within the same region (similar effects, compound heterozygosity), and draws on a measure of “information content” to prioritize results.

Screening for cis-epistatic regions (which may plausibly have evaded selective pressure) prioritizes biologically plausible results while de-emphasizing individual SNPs, which may be significant because of population selection biases, unless they cause exclusively late-onset phenotypes, such as age-related macular degeneration.[86] Avoiding strong model assumptions (additivity, independence) reduces model misspecification biases. Increasing the sample size, instead, does not guard against these biases, so that imposing a higher fixed GWS level in ssGWAS may, somewhat counterintuitively, favor “false positives” over biologically plausible cis-epistatic effects. The main limitation of u-statistics for multivariate data (conceived in the 1940s [87]) is that the amount of memory required became available only with 32-bit operating systems, in 2001, and computations became feasible only with the advent of GPU-enabled cloud computing.

To improve upon the conventional “overly conservative correction”[25] of 7.3, a systematic analysis of GWA studies suggested lowering the GWS level to 7.0 (fixed),[26] and then further by using study-specific empirical approaches.[22] The empirical aGWS decision rule used here accounts for GWAS not being randomized, the absence of a traditional ‘null hypothesis’ in a heritable disease, differences in MAF causing the expected distributions in a QR plot to be convex, and tests in overlapping diplotypes being related.[13]

The combination of a method with higher specificity and a decision strategy with higher sensitivity increased the number of genes above the cut-off while ensuring that the vast majority of genes implicated was related to known pathways in breast cancer etiology, including dysregulation of EEC.

Replication of findings across populations

Conventionally, a lower GWS level required for replication. At the aGWS/2 level, none of most significant ssGWAS results (CGEM: FGFR2, BUB3, MMRN1; EPIC: CHD3, SOHLH2; PBCS: DOCK8) was replicated in another population (Table B of S1 File). Only three genes (AGPAT3, MEGF11, and TRAPPC9) were replicated in both of the other populations, but none for the same SNP. These results are consistent with common lack of replication in ssGWAS.[88] With muGWAS, in contrast, many genes were replicated in at least one population and seven genes were replicated in both of the other populations (FGFR2, TRPM3, AGPAT3, NCOR2, MEGF11, GPC6, and RGS3), although not necessarily in the same intragenic region. Hence, analyses combining the data from several studies (often called “meta-analyses”, even when subject-level data is used) may result in some populations diluting the risk factors present in others.[88]

Our results are consistent with ssGWAS finding recent, highly penetrant mutations, which may differ across populations, while muGWAS has higher power for common cis-epistatic variations, which are more likely to be shared across populations. Even more likely to be shared are genes that carry different variations and different genes with similar contribution to the etiology,[89] consistent with previous findings that breast cancer gene expression signatures have little overlap across populations.[90]

Dysregulated EEC in breast cancer metastasis, angiogenesis, and progression

Genes involved in EEC (e.g., Rab GTPases) are aberrantly expressed in human cancers. [46] Dysregulation of endocytosis-mediated recycling of oncoproteins (e.g., GF receptors and adhesion molecules, including integrins and annexins), can promote progression, migration, and invasion [46, 91]. Cell migration and invasion, which are promoted by EEC of integrins, are also essential features of angiogenesis.[92] In addition, endocytic uptake of lipoproteins is critical for adaptation of cancer to its microenvironment.[93]

Tumor-derived exosomes, 30–150 nm sized extracellular vesicles formed by dysregulated EEC, are critical mediators of intercellular communication between tumor cells and recipient stromal cells in both local and distant microenvironments.[94, 95] Several Rab proteins (Rab2b/5a/9a/27a/27b) are known to function in the selective packaging and production of exosomes in tumor cells (Fig 3, bottom left).[96] Rab27a knockdown in highly metastatic melanoma cells significantly decreased exosome production, primary tumor growth, and metastasis,[97] confirming the role of EEC in generating exosomes.

Dysregulated EEC alters not only exosome biogenesis (vesicular packaging and trafficking), but also the composition of exosomal cargos. Tumor-specific proteins, such as integrins were enriched in exosomes, transferred between cancer cells,[98] and correlated with migration and invasion of recipient cells.[99, 100] Exosome uptake (involving endocytosis [101]) induces non-tumorigenic cells to develop cancer-related phenotypes and the uptake of exosomal integrins promotes migration of these tumor cells.[102]

A recent study revealed that exosomal integrin expression patterns enriched in cancer-derived exosomes involve specific αβ combinations matched to target organs. Proteomic analysis revealed that the exosomal integrin αvβ5 binds to Kupffer cells that mediate liver metastasis, integrins α6β1 and α6β4 are associated with lung metastasis in breast cancer, while integrin β1 (which required for extravasation in metastases [103]) was not organ-specific. [104].

Additionally, other tumor-specific exosomal proteins, such as annexins (calcium-dependent phospholipid-binding proteins known to regulate membrane trafficking and EEC), which are known to correlate with migration and invasion, are also packaged in cancer exosomes [91, 105]. Annexins are frequently overexpressed in breast, liver, prostate, and pancreatic cancers and participate in multiple functions in cancer, including angiogenesis, tumor migration and invasion.[106] In breast cancer, exosomal annexin A2 promotes angiogenesis and vascularization via tissue plasminogen activator (tPA).[91] In pancreatic cancer, exosomal annexin A6 from cancer-associated fibroblasts contributes to tumor cell survival and invasion through annexin A6 / LDL receptor-related protein 1/thrombospondin 1 complex formation.[105]

In summary, EEC plays at least four roles in cancer development; spreading the cancer phenotype horizontally, preparing cancer cells for migration, preparing the distant microenvironment (all via preparation and transmission of exosomes containing integrins), and facilitating migration and invasion (via increasing EEC of integrins). In each case, both endo- and exocytosis are involved, either in donor and target cells or at trailing edge and advancing lamellipodium (Fig 2). Hence, down-regulating “de-railed endocytosis”[46] could have substantial synergistic effects.

The PI cycle in breast cancer

Our findings of PTENP1 (PBCS), TNS1 (EPIC), and SYNJ2 (CGEM) are consistent with known breast cancer mutations in PI3K/PTEN [107, 108] and SYNJ2. That both PI(3,4,5)P3 and PI(3,4)P2 are required to achieve and sustain a malignancy, has been formulated as the “two PI hypothesis”[109] Except for the known PRCKQ, which is regulated by phospholipids via the PI(4,5)P2–PLC–DAG route, however, our analysis identified few genes along the AKT/TSC/mTOR pathway, which is controlled by the “two PIs”. Instead, our results point to EEC, in which virtually all PIs are involved. The closely related set of genes involved in recycling of DAG (DGKQ), influx of PC and PS (ATP8B1, ATP8A1), and influx of LPA and 1-acyl GPI (AGPAT3, AGPAT4) suggests the downregulation of circulating phospholipids as a novel strategy to reduce EEC.

LPA, a known promoter of cell migration and invasion in breast cancer,[110, 111] is produced from LPC by autotaxin (ATX).[112] While ATX mouse knockouts are embryonically lethal, mice that overexpress LPA or ATX develop spontaneous metastatic mammary tumors. A mechanism mediated by G-coupled LPA receptors may cause mesenchymal tumors via endocytosis upregulation involving β-arrestin2 [113] and Arf6.[114]

LPA and LPC in physiologic concentrations have been shown to strongly induce migration of rhabdomyosarcoma (RMS) cells and to be increased by irradiation and chemotherapy in bone marrow.[115] The authors suggested the development of inhibitors of LPA/LPC signaling or “molecules that bind these bioactive lipids” after radio/chemotherapy. However, targeting a single among several redundant receptor/ligand complex may not be sufficiently effective to prevent metastases.[116]

Alkyl-LPCs, which compete with LPC, are in clinical use for treatment of cutaneous metastases in breast cancer, but have shown little activity (and substantial GI side effects) in advanced metastatic breast cancer.[117] From the results presented here, this is consistent with reducing LPC being most effective while cells are still migrating.

As our results suggest, overall EEC upregulation may be caused by multiple variations affecting the PI cycle. Thus, reducing EEC by diminishing the overall phospholipid pool might be a more effective breast cancer treatment than blocking one or even two phosphotransferases, a strategy for which the highly robust PI cycle is designed to compensate. Given the ability of biologic systems to prioritize scarce resources, one would expect this effect to be stronger for tumor cells than for host cells whose functions are routinely prioritized when supplies are scarce. A related approach, substituted myo-inositol (MI) analogues, had already been considered, but was found unlikely to be effective in vivo, because even physiological concentration of MI antagonized the growth inhibitory activity of such analogues.[59]

βCDs are effective in cancer models of migration, invasion, and angiogenesis

A plethora of studies have investigated the effect of methyl-β-cyclodextrin (MβCD) in vitro. For instance, MβCD suppressed translocation of β1 integrin[118] and also invasion activity in three H7 Lewis lung cancer cell lines where highly metastatic cell lines had more β1 integrin.[119]. Breast and prostate cancer cell lines were more sensitive to MβCD-induced cell death than their normal counterparts.[120] In particular, MβCD treatment induced a substantial decrease (40%) in activity of breast cancer resistance protein (BCRP/),[121] which transports PS and PC analogues.[43] In subsequent functional studies, MβCD inhibited spheroid migration and invasion of MDA-MB-241 and ZR751 breast cancer cells [122] and also endocytosis [123] and migration [72] of MCF7 breast cancer cells. MβCD was more toxic for invasive than for non-invasive urothelial cancer cells,[124] and interfered with RTK-[PI2]-PI3K-[PI3]-AKT signaling in HeLa cells.[125] Finally, MβCD reduced breast cancer-induced osteoclast activity in RAW264.7 cells and osteoclastogenic gene expression in MCF-7 cells.[126] Sulfated SβCD also inhibits epithelial cell migration and invasion, but not proliferation [127] and prevents angiogenesis ex vivo in an rat aortic ring assay and an chick embryo collagen onplant assay.[127] The relevance of these in vitro findings was confirmed by several in vivo studies.

MβCD had higher concentration in tumor than in other cells (except kidney and liver involved in its clearance) and reduced tumor volume in mice xenografted with MCF-7 breast cancer or A2780 ovarian carcinoma cells at least as effectively and with less toxicity than doxycyclin,[128] reduced the number of lung metastases in mice implanted with H7-O Lewis lung cancer cells,[119] reduced invasiveness of melanoma,[129] and inhibited growth of primary effusion lymphoma (PEL) in mice.[130] HPβCD was necessary in triple combination treatment for tumor regression in mice implanted with renal cancer cells.[125] and prolonged survival in leukemia mouse models.[131]

βCDs have also seen effective in animal models of several other diseases known to involve endocytosis[132]: Alzheimer’s disease (APP),[133] Parkinson’s disease (α-synuclein),[134] and atherosclerosis (LDL),[135, 136] see Fig 3. However, while HPβCD was well tolerated in most peripheral and central organ systems,[137] it was shown to carry the risk of causing permanent hearing loss in mice,[138] cats,[139, 140] and at least one human.[141] Both intracochlear HPβCD and, in particular, MβCD were seen to be ototoxic in Guinea pigs.[142] This ototoxicity is believed to be due to depriving prestin (SLC26A5) in outer hair cells of cholesterol.[143-145]

Migration and invasion in breast cancer involve cholesterol-unrelated processes

The role of phospholipids emerging from our results, however, suggests a different mechanism than scavenging of cholesterol. This mechanism is consistent with previously reported in vivo results: CAV1 expression in breast cancer stroma increases tumor migration and invasion [146] and CAV1 is required for invadopodia formation specifically by breast cancer cells, where CAV1 knockdown cannot be rescued by cholesterol.[147] Growing MDA-MB-231 breast cancer cells in lipoprotein depleted medium resulted in an 85% decrease in cell migration.[148] LPA activates the Arf6-based mesenchymal pathway for migration and invasion of renal cancer cells, which also originate from cells located within epithelial ductal structures. [114, 143–145]

Limiting the availability PIs would be particularly effective for PI(4)P and PI(4,5)P2 (each at <10%, see Fig 4) and, thus, would likely reduce endocytosis more than lysosomal degradation. In addition, cyclodextrins have been shown to exert their role in NPC treatment by activating rather than downregulating, Ca-dependent lysosomal exocytosis.[149]

From the mechanism of βCD in NPC and elevated cholesterol levels seen in several cancers, including breast cancer,[131] βCDs were thought to reduce cancer growth by lowering cholesterol levels. Early evidence that this might not be the case emerged from the study of exosomes, which play a key role in development of breast cancer.[150, 151] Treatment of MDA-MB-231 breast cancer cells with MβCD inhibited the internalization of exosomes containing integrins,[104] but did so independently of cholesterol.[152]

αCD scavenge phospholipids only, reducing AEs and increasing effectiveness

βCDs is widely believed to act through “cholesterol depletion”,[130, 153] yet βCDs also scavenges phospholipids.[70] From the genetics results, which suggest an overactiv PI cycle (Fig 4) for an age-related decrease of lysosomal function (Fig 3), the effect seen in breast cancer and some of the other diseases may be primarily through scavenging phospholipids. The cavity of αCDs is too small for cholesterol, but large enough for phospholipids.[154, 155] From the in vitro results validating the breast cancer hypothesis generated as part of the U4C challenge (Fig 6), αCDs may be more effective than βCD, yet without the risk of cholesterol-related ototoxicity.

Two types of “controls” have been used: repletion of cholesterol via βCDs “loaded” with cholesterol, and reduction of cholesterol production via statins. Repletion of cholesterol, however, also increases production of phospholipids by freeing acetyl-CoA, the precursor of both phospholipids and cholesterol,[156-158] cholesterol replenishment restores sphingolipid decrease,[118] and statins also lower phospholipids.[159] Hence, neither of these two strategies can “control” against βCDs scavenging phospholipids, rather than cholesterols. Using αCD as a control, however, can answer this question and the above in vitro results suggest that equimolar αCDs are, in fact, at least twice as effective as βCDs, as one would expect if the effect of either CD is caused by its ability to scavenge phospholipids. Hence, our results suggest that many of the previous experiments with βCDs should be redone, this time using αCDs as a control.

αCD is generally recognized as safe (GRAS),(FDA, GRN000155) and approved as an excipient for i.v. alprostadil.[160] Due to higher watersolubility, αCD has lower nephrotoxicity than βCD.[161] HP derivatives of αCD and βCD increase water solubility from 145 and 18.5, respectively, to ≥500 g/L. In mice, the observed ototoxicity order of HPβCD >[p < .002] HPγCD >[p < .02] HPαCD [≈[NS] vehicle] matches the reported order for hemolysis and toxicities in various cell types.[76, 162] In humans, a single dose of up to 3 g/kg/d HPβCD and seven daily doses of 1 g/kg/d were reported to have no adverse effects.[163] In 5-yr old children treated for NPC, 800 mg/kg/d HPβCD i.v. for 12 months was well tolerated.[164]

HPαCD as a potential novel treatment in breast cancer

Given significant redundancy pro-metastatic ligand-receptor complexes, the paradigm of targeting a single receptor-ligand complex has recently been challenged.[116] Although targeting EEC is a promising therapeutic strategy to prevent and treat metastasis,[165] a therapeutic agent is yet to be determined. Our results suggest that metastases in breast cancer rely on upregulation of the highly robust PI cycle and various types of dysregulation along the complex EEC pathway, rather than a simple linear PI pathway. Hence targeting the PI cycle in its entirety may be more effective than targeting individual phosphatases or kinases, or specific genes along the EEC pathway. Cyclodextrins are attractive candidates for a polyvalent approach to treat breast cancer. By modulating several pathways involved in breast cancer, such as altering exosome production and packaging, and impede metastatic colonization, CDs are likely to confer greater protective effects than molecules that have single targets. The selectivity of the smaller αCDs to phospholipids would minimize side effects (e.g., ototoxicity) from βCDs also capturing cholesterol. Given that some CDs are already routinely used clinically, and their pharmacokinetic and toxicity profiles are well established, repeating previous encouraging animal studies of HPβCD, this time using HPαCD could lead rapidly to clinical efficacy trials.

Fig A. Manhattan Plots. Panels are CGEM, PGCS, and EPIC (from top to bottom). Dot size increases from single-SNP (black foreground) to six SNP diplotype (background). Color indicates information content (low: red; high: black). Values that were manually removed, because the diplotype spans LDBs or are higher than overlapping diplotypes are crossed out (white). First line of annotation indicates previous publication implicating this gene; second line implicates genes replicated in at least one of the other two populations.

Fig B: QR-Plot of ssGWAS results by chromosome for CGEM. SNPs too far upstream ("-") or downstream ("+") to be considered related and genes with unknown function (e.g., LOC…, "?") are shown in gray among the results for individual chromosomes and are excluded from the summary plot. The "null" projection (blue) in the summary plot ends at the median among the endpoints of the convex projections for individual chromosomes. Genes to the right of the vertical blue line are above the cut-off for study-specific genome-wide significance.[13].

Fig C: QR-Plot of muGWAS results by chromosome for CGEM. In addition to the annotation used in Fig 2 of S1 File, genes whose significance relies entirely on a single SNP are marked in red ("$") and excluded from the summary plot.

Fig D: String Analysis of Genes aGWS in ssGWAS that are Unrelated to Known Pathways. Connections: Co-Mentioned in PubMed abstracts (green), experimental/biochemical Data (pink), association in curated database (blue);

Table A: Genes Involved in EEC identified in breast cancer. Bold: aGWS. *: from previous GWAS. Underlined: functionally related genes identified in the literature. †: implied.

Table B: Top and replicated genes. Replicated genes (left column) and aGWS results are shown in bold, results below the level of support (row aGWS/2) shown in gray; BC Ref: selected references related to breast cancer; a: Top genes in ssGWAS; b: Top consecutive set of replicated genes by population; c: Other replicated significant genes; d: Additional replicated genes–see Text for details.

1)1-Acylglycerol-3-Phosphate O-Acyltransferase paralogs, convert lysophosphatic acid (LPA) into PA, the second step in de novo phospholipid biosynthesis (hsa00564) 2)Receptor-ligand pair (hsa04060) 3)P-type ATPase paralogs, form flippase complex with TMEM30A, transport amino-phospholipids from the outer to the inner leaflet of membranes 4)homologs, Cowden syndrome, 5)pre/post PA, 6)TRPxn, $quoted in dbGaP from [7].

Table C: Known regulators of integrin cycling in breast cancer. Top: all genes related to breast cancer and endocytosis,[46, 259] *: integrin trafficking genes related to breast cancer.[48].

1)NUMB reduces BC cell migration also by degrading NOTCH;[303] NUMB also interacts with MDM2. MDM2 also induces EMT in breast cancer cells by upregulating Snail.[304].

Knockdown of STXBPx 4 also substantially inhibited β1-integrin recycling in human monocytes.[305].

Fig E: Top and replicated genes.

(PDF)

Click here for additional data file.