  The stellar velocity distribution in the Solar Neighbourhood displays
kinematic substructures, which are possibly signatures of the bar and spiral
arms of the Milky Way and of previous accretion events. These kinematic
substructures -- moving groups -- can be thought of as continuous manifolds in
the 6D phase space, and the ridges in the $V_{\phi}-R$ and arches in the
$V_\phi-V_R$ plane, discovered with the Gaia mission, as projections of these
manifolds. We develop and apply a methodology to perform a blind search for
substructure in the Gaia EDR3 6D data, and obtain a sampling of the manifolds.
The method consists in the execution of the Wavelet Transform in small volumes
of the Milky Way disc, and the grouping of these local solutions into global
ones with a method based on the Breadth-first search algorithm from Graph
Theory. We reveal the complex skeleton of the velocity distribution, sampling
nine main moving groups in a large region of the disc ($6$ kpc, $60$ deg, and
$2$ kpc in the radial, azimuthal, and vertical directions). In the radial
direction, the groups deviate from lines of constant angular momentum that one
would naively expect from first order effect of resonances. The azimuthal
velocity of Acturus, Bobylev, and Hercules is non-axisymmetric. For Hercules,
we measure an azimuthal gradient of $-0.50$ km/s/deg at $R=8$ kpc. We detect a
vertical asymmetry in the azimuthal velocity for Coma Berenices, which is not
expected in a resonance of the bar, supporting previous hypothesis of
incomplete vertical phase-mixing. When we apply the same methodology to
simulations of barred galaxies, we extract substructures corresponding to the
Outer Linbdlad and the 1:1 Resonances and observe patterns consistent with the
data. This data-driven characterization allows for a quantitative comparison
with models, providing a key tool to comprehend the dynamics of the Milky Way.
(Abridged)
