Open access peer-reviewed chapter

Comprehensive Network Analysis of Cancer Stem Cell Signalling through Systematic Integration of Post-Translational Modification Dynamics

Written By

Hiroko Kozuka‐Hata and Masaaki Oyama

Submitted: 11 November 2016 Published: 13 September 2017

DOI: 10.5772/intechopen.69647

From the Edited Volume

Applications of RNA-Seq and Omics Strategies - From Microorganisms to Human Health

Edited by Fabio A. Marchi, Priscila D.R. Cirillo and Elvis C. Mateo

Chapter metrics overview

1,806 Chapter Downloads

View Full Metrics


Post‐translational modifications, such as phosphorylation, acetylation and ubiquitination, are widely known to play various important roles in cellular signalling. Recent significant advances in mass spectrometry‐based proteomics technology enable us not only to comprehensively identify expressed proteins but also to unveil their post‐translational modifications with high sensitivity. In our advanced proteome bioinformatics frameworks, statistical network analyses of large‐scale information on various post‐translational modification dynamics were conducted to define the key machinery for cancer stem cell properties. The bioinformatical approaches using IPA (ingenuity pathway analysis), NetworKIN and a newly developed platform named PTMapper (post‐translational modification mapper) allowed us to perform network‐wide prediction of upstream interactors/kinases with the related information on the diseases and functions, leading to systematic finding of novel drug candidates to regulate aberrant signalling in cancer stem cells. In this chapter, we apply patient‐derived glioblastoma stem cells as a representative model of cancer stem cells to introduce some useful platforms for statistical and mathematical network analyses based on the large‐scale phosphoproteome data.


  • glioblastoma stem cells
  • signal transduction
  • proteomics
  • post‐translational modification
  • network analysis

1. Introduction

Glioblastoma (GBM) is known to be the most common and aggressive brain tumour in adults. Despite the enormous efforts to overcome this tumour for many years, the median survival for GBM patients remains around only 1 year [1]. GBM is characterized by high invasiveness and intratumoral heterogeneity (ITH) [2, 3]. Up to date, it is known that GBM‐ITH contributes to the resistance to chemotherapy, radiation and surgical resection. Since functional diversity is the main feature of multilineage differentiation of cancer stem cells (CSCs) [4, 5], glioblastoma stem cells (GSCs) were thought to be major therapeutic targets of GBM. Furthermore, post‐translational modifications (PTMs) of GSCs are reported to tightly regulate highly tumourigenic potential of GSCs through aberrant signalling [6, 7]. Therefore, it is important to comprehensively elucidate PTM‐based GSC signalling networks for developing the effective treatment of GBM.

Advanced nanoscale liquid chromatography‐tandem mass spectrometry (nanoLC‐MS/MS) enables us to identify and quantify thousands of proteins in a single experiment [8]. Moreover, using the nanoLC‐MS/MS system coupled to the high‐affinity enrichment methods of the peptides with PTMs, we can also acquire in‐depth biological information on PTM dynamics. In this chapter, we introduce high‐resolution shotgun proteomics technology for large‐scale PTM determination in combination with statistical bioinformatics platforms such as IPA [9], NetworKIN [10, 11] and PTMapper [12].


2. System‐wide proteomic analysis of PTM dynamics

PTMs are widely known to play crucial roles in cell fate control, such as proliferation, differentiation and apoptosis. More than 500 kinds of PTMs regarding eukaryotes and prokaryotes have been registered with Unimod, a comprehensive database of protein modifications for mass spectrometry [13]. Recent technological advances in mass spectrometry‐based proteomics in combination with appropriate enrichment techniques for each PTM enable us to perform comprehensive identification and quantification of PTMs [14]. Here, we introduce biochemical purification methods for highly sensitive detection of the representative PTMs: phosphorylation, acetylation and ubiquitination (Figure 1).

Figure 1.

Strategy for mass spectrometry‐based identification of peptides modified with phosphorylation, acetylation and ubiquitination. Regarding ubiquitinated lysine residues, Gly‐Gly remnants are generated from the C‐terminal of ubiquitin as a consequence of tryptic digestion. PTMs: post‐translational modifications, P: phosphorylation, Ac: acetylation, Ub: ubiquitination, TiO2: titanium dioxide.

2.1. Phosphorylation

Protein phosphorylation is recognized as one of the most important and well‐studied PTMs and regulates a variety of biological processes by transmitting diverse external signals [15, 16]. About as many as 280,000 phosphorylation sites have already been registered in PhosphoSitePlus, a knowledgebase containing non‐redundant mammalian PTMs [17]. Titanium dioxide (TiO2), which has very high affinity for phosphorylated peptides, is widely used for large‐scale phosphoproteome analysis [18, 19].

2.2. Acetylation

Lysine acetylation plays a key role in modulating transcriptional regulation through the coordinated function of histone acetyltransferases (HATs) and histone deacetylases (HDACs) [20]. The stabilization of p53, one of the most important transcription factors, is reported to greatly depend on lysine acetylation [21]. Thousands of lysine acetylation sites can be identified using an antibody against acetyl‐lysine in combination with a high‐resolution mass spectrometry system [22, 23].

2.3. Ubiquitination

The ubiquitin system transmits protein degradation signal to proteasome as well as regulates multiple cellular functions such as cell‐cycle progression, DNA repair and transcriptional regulation. Dysfunction of this system leads to various pathological conditions [24]. Ubiquitination sites are detected as diglycine (Gly‐Gly) remnants on the modified lysine residues, which are generated by tryptic digestion of ubiquitinated proteins [25, 26].


3. Systematic characterization of the phosphoproteome dynamics in GSCs

The quantitative information on the phosphoproteome dynamics can provide us with systematic description of the key machinery for cellular signalling. In this section, we introduce two examples of global phosphoproteome analyses of GSCs using SILAC (stable isotope labelling by amino acids in cell culture)‐based quantitative technique [27, 28] (Figure 2). One was carried out using epidermal growth factor (EGF) to elucidate the mechanism for stemness maintenance of GSCs [29], whereas the other was conducted through serum‐induced differentiation of GSCs to unveil the key pathways responsible for disrupting stemness characteristics [30].

Figure 2.

Schematic workflow for quantitative proteome analysis using SILAC, a representative relative quantitation technique based on metabolic labelling of specific amino acids such as arginine. Two populations of GSCs were cultured in the media supplemented with 12C614N4‐Arg (light) or 13C615N4‐Arg (heavy), respectively. After one of the two cell populations was stimulated/perturbed, both of the cells were lysed, equally combined and enzymatically digested to perform nanoLC‐MS/MS analyses. The intensity of each mass peak is used for relative quantitation of each peptide with high accuracy.

3.1. Global quantitative phosphoproteome analyses of EGF‐stimulated GSCs

EGF is known to be essential for maintenance and growth of GSCs [31]. The quantitative phosphoproteomic analysis of EGF‐stimulated GSCs was performed to acquire network‐wide information on the molecules related to stemness maintenance. As a result, a total of 6073 phosphopeptides from 2282 phosphorylated proteins were identified, leading to quantitative classification of 516 upregulated and 275 downregulated phosphorylation sites [29].

3.1.1. IPA‐based network analysis

IPA canonical pathway analysis was then performed using SILAC‐based quantitative phosphoproteome data on EGF‐stimulated GSCs [29] (Figure 3). Protein synthesis‐related pathways (EIF2 signalling, mTOR signalling) and cell cycle regulation‐related pathways (cyclins and cell cycle regulation, cell cycle: G1/S checkpoint regulation, cell cycle: G2/M DNA damage checkpoint regulation) were extracted with statistical significance (‐log (p‐value) > 5).

Figure 3.

IPA‐based pathway analysis of the quantitative phosphoproteome data on EGF‐stimulated GSCs. (A) The significant canonical pathways across the entire dataset (‐log (p‐value) > 5). (B) The mTOR signalling pathway is representatively depicted with the predicted information on the biological activities related to this pathway.

3.1.2. Upstream kinase prediction analysis

Protein phosphorylation is known to be controlled by specific kinases depending on consensus sequence motifs of substrates [32]. The motif‐x algorithm [33, 34] is applicable to statistical extraction of significant consensus sequence motifs from the large‐scale phosphoproteome data on EGF‐stimulated GSCs (Figure 4(A) and (B)).

Figure 4.

Phosphorylation site‐oriented network analysis of the quantitative phosphoproteome data on EGF‐stimulated GSCs. The consensus sequence motifs surrounding the quantitatively regulated phosphorylation sites regarding (A) downregulation and (B) upregulation can be described as a result of the motif‐x analyses. (C) The numerical distribution of the putative kinases predicted by NetworKIN. The colour of cells reflects the number of the predicted kinases for each consensus sequence as described in (A) and (B).

NetworKIN [10, 11] is designed to predict upstream kinases based on the sequence motifs around the functionally regulated phosphorylation sites through construction of the related protein‐protein interaction (PPI) networks using STRING [35]. The NetworKIN algorithm enables further interpretation of the results obtained from the motif‐x analyses (Figure 4 (C)).

3.2. Global quantitative phosphoproteome analyses of serum‐induced GSCs

CSCs are regarded as one of the most clinically important cell populations in causing tumour heterogeneity, which is responsible for the resistance to chemotherapy [36]. As recent studies have demonstrated that non‐CSCs can also readily acquire CSC‐like characteristics [37], it is very important to figure out the detailed mechanisms underlying CSC differentiation and understand the principle of their heterogeneity. Serum‐induced phosphoproteome dynamics in GSCs was measured to systematically elucidate the regulatory nodes for stemness alteration over the entire signalling networks [30]. Among 2876 phosphorylation sites on 1584 proteins identified, 732 phosphorylation sites on 419 proteins were found to be regulated through serum‐induced differentiation. The integrative network analyses of the quantitative phosphoproteome data using various bioinformatical tools including IPA and NetworKIN indicated that transforming growth factor‐β receptor type‐2 (TGFBR2) might be one of the crucial upstream regulators concerning GSC alteration (Figure 5).

Figure 5.

Upstream kinase/regulator analyses based on the regulated phosphoproteome data on serum‐induced GSCs. (A) Heatmap of the over‐representation p‐values calculated for each predicted kinase using PhosphoSiteAnalyzer, a bioinformatical platform for the NetworKIN prediction results from the phosphoproteome data [38]. The subset ‘serum (−)’ indicates SILAC ratio > 2.0, whereas ‘serum (+)’ shows SILAC ratio < 0.5. TGFBR2 and ACVR2A/B‐specific phosphorylation sites were predicted to be significantly enriched in the ‘serum (−)’ subset (adjusted p‐value < 0.05). (B) Upstream regulator analysis by IPA. The top 10 upstream regulators relevant to the regulated phosphoproteome are shown with the corresponding score (−log [p‐value]). (C) IPA‐based description of TGF‐β1 and the target molecules in the phosphoproteome data. Dashed lines represent indirect interactions caused by TGF‐β1, adapted from Ref. [30].


4. Development of advanced bioinformatical platforms for complicated kinase‐substrate interaction networks

Although shotgun proteomics strategy based on advanced nanoLC‐MS/MS system can provide us with large‐scale information on various kinds of PTMs, there are only a few PTM‐based network analysis tools available compared to conventional protein‐protein interaction (PPI). Recently, CEASAR: connecting enzymes and substrates at amino acid resolution [39] and PhosphoPath [40] were developed to visualize kinase‐substrate interactions in a phosphorylation site‐oriented manner. CEASAR was designed to provide a high‐resolution map of kinase‐phosphorylation networks based on functional protein microarrays and bioinformatics analysis. On the other hand, PhosphoPath was developed as a Cytoscape app [41] to visualize both quantitative proteome and phosphoproteome data using PPI information extracted from BioGRID [42] and PhosphoSitePlus [17]. Recently, we also have developed a Cytoscape‐based bioinformatical platform named ‘post‐translational modification mapper (PTMapper)’ to visualize kinase‐substrate interactions regarding multiple phosphorylation sites on signalling molecules (Figure 6) [12]. The kinase‐phosphorylation site interaction dataset for this platform was integratively generated from PhosphoSitePlus [17], Phospho.ELM [43], PhosphoNetworks [44] and Uniprot KB [45], leading to construction of phosphorylation site‐oriented PPI networks using Pathway Commons [46]. We applied this platform to extract crucial kinase‐substrate interactions from the quantitative phosphoproteome data on EGF‐stimulated GSCs [29]. As a result, p70S6K and Lyn were significantly extracted as key regulators (Figure 7).

Figure 6.

Construction of phosphorylation‐oriented PPI networks via PTMapper. (A) Workflow for the visualization of kinase‐phosphorylation site relationships in PPI networks via PTMapper. Phosphorylation sites are connected with the parental protein nodes in PPI networks and the upstream kinases are then added to the phosphorylation sites. (B) Phosphorylation site-oriented networks constructed from the phosphoproteome data on EGF-stimulated glioblastoma stem cells. The solid arrows represent functionally directed protein‐protein interactions or kinase‐substrate interactions, whereas the dotted lines show the linkages of proteins and their phosphorylation sites, adapted from Ref. [12].

Figure 7.

Comparison of the sub‐networks extracted from EGF‐dependent phosphorylation dynamics of glioblastoma stem cells. (A) Schematic procedure for the evaluation of PTMapper‐based network construction. (B) The most significantly regulated sub‐networks extracted from the conventional protein interaction network. (C) The phosphorylation site-oriented network generated via PTMapper. The nodes surrounded by the border with the upper-right numbers indicate the common molecules in the two types of the sub-networks. The solid arrows represent functionally directed protein-protein interactions or kinase-substrate interactions, whereas the dotted lines show the linkages of proteins and their phosphorylation sites. The dashed circles indicate p70S6K and Lyn, adapted from Ref. [12].


5. Perspectives and conclusions

The bioinformatical description of GSC signalling dynamics based on the global quantitative phosphoproteome data led to network‐wide extraction of critical molecules and their related pathways for defining stemness characteristics. Further integrative description of multiple PTM dynamics in GSCs will deepen our understanding of the nature of their cell signalling complexity at the network level. We believe that shotgun proteomics‐based quantitative analyses of cancer stem cell signalling networks in combination with various statistical and mathematical platforms will pave the way to establish new directions towards systematic evaluation of drug targets in a cell‐type specific manner.



We thank Dr. Yuta Narushima for his technical support. We are also thankful to all the members of Medical Proteomics Laboratory, The Institute of Medical Science, The University of Tokyo. This work was supported by grants‐in‐aid for scientific research on innovative areas (integrative understanding of biological signalling networks based on mathematical science) and grant‐in‐aid for scientific research (C).


  1. 1. Furnari FB, et al. Malignant astrocytic glioma: Genetics, biology, and paths to treatment. Genes & Development. 2007;21(21):2683-2710. DOI: 10.1101/gad.1596707
  2. 2. Kreso A, Dick JE. Evolution of the cancer stem cell model. Cell Stem Cell. 2014;14(3):275-291. DOI: 10.1016/j.stem.2014.02.006
  3. 3. Burrell RA, McGranahan N, Bartek J, Swanton C. The causes and consequences of genetic heterogeneity in cancer evolution. Nature. 2013;501(7467):338-345. DOI: 10.1038/nature12625
  4. 4. Vescovi AL, Galli R, Reynolds BA. Brain tumour stem cells. Nature Reviews Cancer. 2006;6(6):425-436. DOI: 10.1038/nrc1889
  5. 5. Qazi MA, Vora P, Venugopal C, Sidhu SS, Moffat J, Swanton C, Singh SK. Intratumoral heterogeneity: Pathways to treatment resistance and relapse in human glioblastoma. Annals of Oncology. 2017. DOI: 10.1093/annonc/mdx169
  6. 6. Takebe N, Harris PJ, Warren RQ, Ivy SP. Targeting cancer stem cells by inhibiting Wnt, Notch, and Hedgehog pathways. Nature Reviews Clinical Oncology. 2011;8(2):97-106. DOI: 10.1038/nrclinonc.2010.196
  7. 7. Wurdak H, et al. An RNAi screen identifies TRRAP as a regulator of brain tumor‐initiating cell differentiation. Cell Stem Cell. 2010;6(1):37-47. DOI: 10.1016/j.stem.2009.11.002
  8. 8. Aebersold R, Mann M. Mass spectrometry‐based proteomics. Nature. 2003;422(6928):198-207. DOI: 10.1038/nature01511
  9. 9. Krämer A, Green J, Pollard Jr J, Tugendreich S. Causal analysis approaches in ingenuity pathway analysis. Bioinformatics. 2014;30(4):523-530. DOI: 10.1093/bioinformatics/btt703. Available from:
  10. 10. Horn H, et al. KinomeXplorer: An integrated platform for kinome biology studies. Nature Methods. 2014;11(6):603-604. DOI: 10.1038/nmeth.2968. Available from:
  11. 11. Linding R, et al. Systematic discovery of in vivo phosphorylation networks. Cell. 2007;129(7):1415-1426. DOI: 10.1016/j.cell.2007.05.052
  12. 12. Narushima Y, Kozuka‐Hata H, Tsumoto K, Inoue J, Oyama M. Quantitative phosphoproteomics‐based molecular network description for high‐resolution kinase‐substrate interactome analysis. Bioinformatics. 2016;32(14):2083-2088. DOI: 10.1093/bioinformatics/btw164. Available from:‐narushima/PTMapper/
  13. 13. Creasy DM, Cottrell JS. Unimod: Protein modifications for mass spectrometry. Proteomics. 2004;4(6):1534-1536. DOI: 10.1002/pmic.200300744. Available from:
  14. 14. Mann M, Jensen ON. Proteomic analysis of post‐translational modifications. Nature Biotechnology. 2003;21(3):255-261. DOI: 10.1038/nbt0303‐255
  15. 15. Hunter T. Protein kinases and phosphatases: The yin and yang of protein phosphorylation and signalling. Cell. 1995;80(2):225-236
  16. 16. Hunter T. Signalling‐‐2000 and beyond. Cell. 2000;100(1):113-127
  17. 17. Hornbeck PV, et al. PhosphoSitePlus, 2014: Mutations, PTMs and recalibrations. Nucleic Acids Research. 2015;43(Database issue):D512‐D520. DOI: 10.1093/nar/gku1267. Available from:
  18. 18. Larsen MR, Thingholm TE, Jensen ON, Roepstorff P, Jørgensen TJ. Highly selective enrichment of phosphorylated peptides from peptide mixtures using titanium dioxide microcolumns. Molecular & Cellular Proteomics. 2005;4(7):873-886. DOI: 10.1074/mcp.T500007‐MCP200
  19. 19. Olsen JV, et al. Global, in vivo, and site‐specific phosphorylation dynamics in signaling networks. Cell. 2006;127(3):635-648. DOI: 10.1016/j.cell.2006.09.026
  20. 20. Kouzarides T. Acetylation: A regulatory modification to rival phosphorylation? EMBO Journal. 2000;19(6):1176-1179. DOI: 10.1093/emboj/19.6.1176
  21. 21. Yang XJ, Seto E. Lysine acetylation: Codified crosstalk with other posttranslational modifications. Molecular Cell. 2008;31(4):449-461. DOI: 10.1016/j.molcel.2008.07.002
  22. 22. Kim SC, et al. Substrate and functional diversity of lysine acetylation revealed by a proteomics survey. Molecular Cell. 2006;23(4):607-618. DOI: 10.1016/j.molcel.2006.06.026
  23. 23. Choudhary C, et al. Lysine acetylation targets protein complexes and co‐regulates major cellular functions. Science. 2009;325(5942):834-840. DOI: 10.1126/science.1175371
  24. 24. Hershko A, Ciechanover A. The ubiquitin system. Annual Review of Biochemistry. 1998;67:425-479. DOI: 10.1146/annurev.biochem.67.1.425
  25. 25. Xu G, Paige JS, Jaffrey SR. Global analysis of lysine ubiquitination by ubiquitin remnant immunoaffinity profiling. Nature Biotechnology. 2010;28(8):868-873. DOI: 10.1038/nbt.1654
  26. 26. Wagner SA, et al. A proteome‐wide, quantitative survey of in vivo ubiquitylation sites reveals widespread regulatory roles. Molecular & Cellular Proteomics. 2011;10(10):M111.013284. DOI: 10.1074/mcp.M111.013284
  27. 27. Ong SE, et al. Stable isotope labelling by amino acids in cell culture, SILAC, as a simple and accurate approach to expression proteomics. Molecular & Cellular Proteomics. 2002;1(5):376-386
  28. 28. Blagoev B, et al. A proteomics strategy to elucidate functional protein‐protein interactions applied to EGF signaling. Nature Biotechnology. 2003;21(3):315-318. DOI: 10.1038/nbt790
  29. 29. Kozuka‐Hata, et al. Phosphoproteome of human glioblastoma initiating cells reveals novel signaling regulators encoded by the transcriptome. PLoS One. 2012;7(8):e43398. DOI: 10.1371/journal.pone.0043398
  30. 30. Narushima Y, et al. Integrative network analysis combined with quantitative phosphoproteomics reveals transforming growth Factor‐beta receptor type‐2 (TGFBR2) as a novel regulator of glioblastoma stem cell properties. Molecular & Cellular Proteomics. 2016;15(3):1017-1031. DOI: 10.1074/mcp.M115.049999
  31. 31. Pollard SM, et al. Glioma stem cell lines expanded in adherent culture have tumor‐specific phenotypes and are suitable for chemical and genetic screens. Cell Stem Cell. 2009;4(6):568-580. DOI: 10.1016/j.stem.2009.03.014
  32. 32. Seet BT, Dikic I, Zhou MM, Pawson T. Reading protein modifications with interaction domains. Nature Reviews Molecular Cell Biology. 2006;7(7):473-483. DOI: 10.1038/nrm1960
  33. 33. Schwartz D, Gygi SP. An iterative statistical approach to the identification of protein phosphorylation motifs from large‐scale data sets. Nature Biotechnology. 2005;23(11):1391-1398. DOI: 10.1038/nbt1146. Available from: http://motif‐
  34. 34. Chou MF, Schwartz D. Biological sequence motif discovery using motif‐x. Current Protocols in Bioinformatics. 2011;13(13):15-24. DOI: 10.1002/0471250953.bi1315s35
  35. 35. Szklarczyk D, et al. The STRING database in 2017: Quality‐controlled protein‐protein association networks, made broadly accessible. Nucleic Acids Research. 2017;45(D1):D362‐D368. DOI: 10.1093/nar/gkw937. Available from: http://string‐
  36. 36. Meacham CE, Morrison SJ. Tumour heterogeneity and cancer cell plasticity. Nature. 2013;501(7467):328-337. DOI: 10.1038/nature12624
  37. 37. Chaffer CL, et al. Poised chromatin at the ZEB1 promoter enables breast cancer cell plasticity and enhances tumorigenicity. Cell. 2013;154(1):61-74. DOI: 10.1016/j.cell.2013.06.005
  38. 38. Bennetzen MV, Cox J, Mann M, Andersen JS. PhosphoSiteAnalyzer: A bioinformatic platform for deciphering phosphoproteomes using kinase predictions retrieved from NetworKIN. Journal of Proteome Research. 2012;11(6):3480-3486. DOI: 10.1021/pr300016e. Available from:
  39. 39. Newman RH, et al. Construction of human activity‐based phosphorylation networks. Molecular Systems Biology. 2013;9:655. DOI: 10.1038/msb.2013.12
  40. 40. Raaijmakers LM, et al. PhosphoPath: Visualization of Phosphosite‐centric dynamics in temporal molecular networks. Journal of Proteome Research. 2015;14(10):4332-4341. DOI: 10.1021/acs.jproteome.5b00529. Available from:
  41. 41. Shannon P, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Research. 2003;13(11):2498-2504. DOI: 10.1101/gr.1239303. Available from:
  42. 42. Chatr‐Aryamontri A, et al. The BioGRID interaction database: 2015 update. Nucleic Acids Research. 2015;43(Database issue):D470‐D478. DOI: 10.1093/nar/gku1204. Available from:
  43. 43. Dinkel H, et al. Phospho.ELM: A database of phosphorylation sites‐‐update 2011. Nucleic Acids Research. 2011;39(Database issue):D261‐D267. DOI: 10.1093/nar/gkq1104. Available from:
  44. 44. Hu J, et al. PhosphoNetworks: A database for human phosphorylation networks. Bioinformatics. 2014;30(1):141-142. DOI: 10.1093/bioinformatics/btt627. Available from:
  45. 45. Magrane M, et al. UniProt Knowledgebase: A hub of integrated protein data. Database (Oxford). 2011;bar009. DOI: 10.1093/database/bar009. Available from:
  46. 46. Cerami EG, et al. Pathway commons, a web resource for biological pathway data. Nucleic Acids Research. 2011;39(Database issue):D685‐D690. DOI: 10.1093/nar/gkq1039. Available from:

Written By

Hiroko Kozuka‐Hata and Masaaki Oyama

Submitted: 11 November 2016 Published: 13 September 2017