How does one prove that the lowest-order nonvanishing multipole moment of a charge distribution is independent of the origin, for arbitrary $\ell$?

The multipole moments of a distribution are independent of origin if all the lower terms are zero. I can explicitly verify this statement by hand up to the quadrupole level, but is there any straight forward way of proving this statement?

Generally speaking, whatever definition of a mulitpole moment you use, it will look something like this: $$Q_m^{(\ell)} = \int x_{i_1}x_{i_2}\cdots x_{i_\ell} \,c^{i_1,\cdots,i_\ell}_{m} \, \rho(\mathbf r) \mathrm d\mathbf r = \int p_m^{(\ell)}(\mathbf r) \, \rho(\mathbf r) \mathrm d\mathbf r$$ i.e. as the charge-weighted integral of some homogeneous polynomial $$p_m^{(\ell)}(\mathbf r) = c^{i_1,\cdots,i_\ell}_{m} \, x_{i_1}x_{i_2}\cdots x_{i_\ell}$$ of degree $\ell$ with coefficients $c^{i_1,\cdots,i_\ell}_{m}$ (with Einstein summation understood). If you displace the origin on such a homogeneous polynomial, then you will get a messier polynomial for which the leading-order terms are unchanged but which now has a bunch of lower-order polynomials contributing: $$p_m^{(\ell)}(\mathbf r+\mathbf r_0) = p_m^{(\ell)}(\mathbf r) + \sum_{k=0}^{\ell-1}f^{(\ell-k)}(\mathbf r)g^{(k)}(\mathbf r_0),$$ where $f^{(\ell-k)}(\mathbf r)$ and $g^{(k)}(\mathbf r_0)$ are homogeneous polynomials of degree $\ell-k$ and $k$ respectively, which then transforms the multipole moments as \begin{align} Q_m^{(\ell)} \mapsto \tilde Q_m^{(\ell)} & = \int p_m^{(\ell)}(\mathbf r+\mathbf r_0) \, \rho(\mathbf r) \mathrm d\mathbf r \\ & = \int p_m^{(\ell)}(\mathbf r) \, \rho(\mathbf r) \mathrm d\mathbf r + \sum_{k=0}^{\ell-1}g^{(k)}(\mathbf r_0) \int f^{(\ell-k)}(\mathbf r)\, \rho(\mathbf r) \mathrm d\mathbf r. \end{align}
And now here is the rub: the lower-order polynomials $f^{(\ell-k)}(\mathbf r)$ are all linear combinations of the $p_{m'}^{(\ell-k)}(\mathbf r)$ that form the integral kernels for the lower-order multipole moments, which means that they must vanish when integrated.
The problem, of course, is to actually show that the $f^{(\ell-k)}(\mathbf r)$ are linear combinations of the $p_{m'}^{(\ell-k)}(\mathbf r)$, because the latter are not full basis of the homogeneous polynomial space: as a simple example, if you look at the quadrupole layer, the $p_{m'}^{(\ell-k)}(\mathbf r)$ read basically as $$\{xy,xz,yz, x^2-y^2, 2z^2-x^2-y^2\},$$ with five members instead of six ─ they are missing the homogeneous polynomial $x^2+y^2+z^2$ from their span.
That, then, gives the real question of the proof: do the $f^{(\ell-k)}(\mathbf r)$ fall inside the span of this subspace? And this is where the nebulous edges of our original question come back to bite us, because if we started off using the wrong (i.e. non-traceless) definition of the multipole moments, then those encode too much information (i.e. more than actually gets used in the multipole expansion of the field) and those can end up producing $f^{(\ell-k)}(\mathbf r)$s that fall outside of the relevant subspace; for those choices of multipole moment, the moments themselves may not be origin independent. (As an example, consider the non-traceless octupole moment $\int x^3 \rho(\mathbf r)\mathrm d\mathbf r$.)
So how do you fix that? It comes down to the traceless-ness-ing procedure, which isn't treated in full depth for arbitrary order $\ell$ in the references I can get my hands on at the moment, so I'll leave it to you to formalize once you have a strict definition of the multipole moments that you want to use. But hopefully this answer will make it clearer how the roadmap goes.