The shock response of CuxZr100-x (x = 30, 50 and 70) metallic glasses (MGs) is understood using large-scale molecular dynamics simulations. Piston velocities from Up = 0.125 to 2.5 km/s are simulated, corresponding to shock pressures from 3 to 130 GPa. For each composition, the metallic glasses exhibit different shock wave propagation regimes: (1) single elastic shock wave for Up < 0.25 km/s, (2) split elastic and plastic shock waves for 0.25 < Up 0.75 km/s. Within the split wave and overdriven regimes, the amplitude of the elastic precursor increases with increasing shock intensity, thereby indicating a pressure-dependent yield criterion. Hugoniot states are strongly dependent on the Cu content of the MG with Cu70Zr30 exhibiting higher resistance to plastic deformation than either Cu50Zr50 or Cu30Zr70. Plastic deformation in the MGs initiates via formation of localized regions of high von Mises shear strain, known as shear transformation zones (STZs). At low impact velocities, but above the Hugoniot elastic limit, STZ nucleation is disperse behind the shock front. As impact velocity is increased, STZ nucleation becomes more homogeneous, eventually leading to shock-induced melting, which is identified in this work via high atomic diffusivity. Both the flow stress in the plastically deformed MG and the critical shock pressure associated with melting increase with increasing Cu content over the entire range of piston velocities. The evolution of short-range order in the MG samples is analysed using a polydisperse Voronoi tessellation method. Cu-centered polyhedra with full icosahedral symmetry are most resistant to evolution under shock loading, independent of MG composition. The collapse of nanosized voids within the MG during shock is also examined.