We present the first three-dimensional (3D) simulations of the large-scale mixing that takes place in the shock-heated stellar layers ejected in the explosion of a 15.5 solar-mass blue supergiant star. The outgoing supernova shock is followed from its launch by neutrino heating until it breaks out from the stellar surface more than two hours after the core collapse. Violent convective overturn in the post-shock layer causes the explosion to start with significant asphericity, which triggers the growth of Rayleigh-Taylor (RT) instabilities at the composition interfaces of the exploding star. Deep inward mixing of hydrogen (H) is found as well as fast-moving, metal-rich clumps penetrating with high velocities far into the H-envelope of the star as observed, e.g., in the case of SN 1987A. Also individual clumps containing a sizeable fraction of the ejected iron-group elements (up to several 0.001 solar masses) are obtained in some models. The metal core of the progenitor is partially turned over with Ni-dominated fingers overtaking oxygen-rich bullets and both Ni and O moving well ahead of the material from the carbon layer. Comparing with corresponding 2D (axially symmetric) calculations, we determine the growth of the RT fingers to be faster, the deceleration of the dense metal-carrying clumps in the He and H layers to be reduced, the asymptotic clump velocities in the H-shell to be higher (up to ~4500 km/s for the considered progenitor and an explosion energy of 10^{51} ergs, instead of <2000 km/s in 2D), and the outward radial mixing of heavy elements and inward mixing of hydrogen to be more efficient in 3D than in 2D. We present a simple argument that explains these results as a consequence of the different action of drag forces on moving objects in the two geometries. (abridged)