The galactocentric radial metallicity gradients of GC systems can be used to probe the assembly histories of galaxies but doing so is complicated by: 1.) a GC system is an aggregate population of several (or many) different subpopulations of different origins, and 2.) computing the gradients with respect to projected galactocentric distances can artificially flatten the gradient. I describe a hierarchical statistical framework measure the radial metallically gradients of the GC subpopulations and its application the brightest cluster galaxy M87 using spectroscopically-derived metallicities for a large sample of GCs. Contrary to much previous work I find that the metal-rich and metal-poor subpopulations have different gradients, indicating they were affected by separate events. Moreover, I find the gradients change between the inner and outer halo, indicating that M87 has at least two merger events in its assembly history. I describe the results of doing ``weak chemical tagging'' on physically-motivated spectral stacks of the GC subpopulations and compare these results to the stellar parameters derived spatially-resolved spectroscopy for the galaxy light of M87. I will discuss these results in context of the formation and assembly history of M87 and assess how well predictions from cosmological galaxy formation simulations match our results.