This same question has been asked in the physics section and answered with the usual physics proof, which is unsatisfying to me as I wish to derive the symmetry properties DIRECTLY FROM THOSE OF THE DIFFERENTIAL OPERATOR. Here are two of these proofs that hide the true reason for invariance.
The equation at stake (I would gladly take a general differential operator $L$ but let us start simple) $$\text{Klein-Gordon}\qquad \left\lbrace \begin{aligned} \big(\square + m^2\big)\, \Delta(x,y) &= 0\\ \Delta\left( x^0, \mathbf{x}, x^0, \mathbf{y} \right) & = 0\\ \lim_{t\to x_0} \partial_t\, \Delta\left( t, \mathbf{x}, x^0,\mathbf{y}\right) & = -i \hbar\, \delta^{(3)}(\mathbf{x}-\mathbf{y}) \end{aligned}\right. \label{1}\tag{1} $$ Remark: $\Delta(x,y)$ is a solution (for each of its variable $x$ or $y$) of the homogeneous equation. I call it a propagator (namely Pauli-Jordan or commutator function) rather than a Green function which I leave to designate a solution of the equation with a $\delta(x-y)$ on the r.h.s.. The initial conditions are such that $$ i \hbar\, \Delta(x,y) := \left\langle 0\! \mid \left[ \varphi(x), \varphi(y) \right] \mid\! 0 \right\rangle \qquad \text{where}\quad \left[ \varphi(x), \varphi(y) \right] := \varphi(x)\cdot \varphi(y)- \varphi(y)\cdot \varphi(x)$$ where $\varphi$ are operator-valued solutions of (\ref{1}) satisfying precisely the initial conditions at "equal-time".
These operators act on a Hilbert space $\mathcal{H}$ which admits a unitary representation $U:\mathcal{P} \to \mathcal{U}(\mathcal{H})$ of the Poincaré group. Moreover $\mathcal{H}$ is postulated to contain an invariant vector $\mid 0\rangle$ (a.k.a. vacuum state): $ \displaystyle\forall\ (w,\Lambda)\in\mathcal{P},\; U(w,\Lambda) \mid\! 0\rangle =\, \mid\! 0\rangle$. "Dually", $\mathcal{P}$ acts on the field by conjugation: $\displaystyle U(w,\Lambda)^{-1}\!\cdot\hspace{.5pt} \varphi(x)\cdot U(w,\Lambda) = \varphi\left( \Lambda\cdot x + w\right)$. One has $$ \begin{split} \Delta\left( \Lambda\cdot x + w,\hspace{.7pt} \Lambda\cdot y + w\right) &= \left\langle 0 \left\lvert\, \left[ U(w,\Lambda)^{-1}\!\cdot\hspace{.5pt} \varphi(y)\cdot U(w,\Lambda) ,\hspace{.7pt} U(w,\Lambda)^{-1}\!\cdot\hspace{.5pt} \varphi(x)\cdot U(w,\Lambda) \right]\, \right\rvert 0 \right\rangle\\ &= \left\langle 0 \left\lvert\ U(w,\Lambda)^{-1}\!\cdot\hspace{.5pt} \left[ \varphi(x) , \hspace{.5pt} \varphi(y) \right] \cdot U(w,\Lambda) \, \right\rvert 0 \right\rangle = \Delta (x,y) \end{split} \label{2}\tag{2} $$
The previous proof uses (in fact, postulates) many things, but most importantly that the "quantized" field satisfy the same equation as the "classical" one and could probably be brought to the following one, namely that one can derive an explicit expression for $ \Delta(x,y)$: $\varphi$ is a solution of the homogeneous equation. Its Fourier transform $\hat{\varphi}$ thus satisfies (cf. also here)
$$ \left(-k^2 +m^2\right) \hat{\varphi}(k)= 0 \qquad\text{where}\quad k^2:= (k^0)^2 - \left\lVert \mathcal{k}\right\rVert^2 \label{3}\tag{3}$$
which imposes that $\hat{\varphi}$ vanishes as soon as $\left(-k^2 +m^2\right)\neq 0$. It thus has support in the zero locus of this polynomial which corresponds to a hyperboloic, and is usually parametrized as the graph of the two functions $$ \left\lbrace \begin{aligned} \mathbb{R}^3 & \longrightarrow \enspace \mathbb{R}\\ \mathbf{k}\ & \longmapsto \pm\, \omega_{\mathbf{k}} := \pm \sqrt{\left\lVert \mathcal{k}\right\rVert^2 + m^2} \end{aligned} \right. \label{4}\tag{4}$$
All in all, one conventionnally decomposes the Klein-Gordon field as $$ \varphi(x) = \int_{\mathbb{R}^3} \left(\rule{0mm}{10pt} a_{\mathbf{k}}\, e^{-i k\cdot x} + a^{\dagger}_{\mathbf{k}}\, e^{i k\cdot x} \right) \frac{ d^3 \mathbf{k}}{(2\pi)^3 \sqrt{2\hspace{.4mm} \omega_{\mathbf{k}}}} \quad\text{where} \quad \left\lbrace \begin{aligned} k\cdot x &:= k^0\, x^0 - \mathbf{k}\cdot \mathbf{x}\\ k^0 &:= \omega_{\mathbf{k}}\end{aligned} \right.\label{5}\tag{5}$$
Remark: The point of this decomposition is that the space-time evolution has been factorized in the exponential terms. There are several conventions of this decomposition...
Going back to the quantized case (which though it manipulates operator, has the same space-time dependence and which in the end is brought to a real valued distribution by evaluating with a linear form (state)) by using the commutation relation for the so-called annihilation and creation
$$ \left[ \mathbf{a}_{\mathbf{h}},\mathbf{a}^{\dagger}_{\mathbf{k}}\right] =(2\hspace{.7pt}\pi)^3\, \delta^{(3)}(\mathbf{h}-\mathbf{k}) \label{6}\tag{6}$$
One finally obtains (cf. "An introduction to QFT" Peskin and Schroeder (2.29) p.21, (2.50) & (2.53) p.27-28)
$$i \, \Delta(x,y) = \int_{\mathbb{R}^3} \left( e^{-i k\cdot (x-y)}- e^{ik\cdot (x-y)} \right) \frac{d^3 \mathbf{k}}{(2\pi)^3\, 2\hspace{.9pt} \omega_{\mathbf{k}}} \label{7}\tag{7}$$
which is invariant when one simultaneously translated $x$ and $y$ (because the dependence in $x,y$ is of the form $x-y$) and when one simultaneously apply a Lorentz transform as the measure $\displaystyle \frac{d^3 \mathbf{k}}{\omega_{\mathbf{k}}}$ is Lorentz invariant and that one integrales over the whole (positive energy) hyperboloid
What I wish is a simple/direct/elementary/general proof going along the following line: let $A\in \mathrm{GL}_n(\mathbb{R}),\ M:\left\lbrace \begin{aligned} G & \longrightarrow \mathrm{GL}_n(\mathbb{R})\\ g\, & \longmapsto\ M_g \end{aligned}\right.$ a representation of a group $G$ such that $\forall\ g\in G,\; [A,M_g]=0$. Then one has for the inverse $B:= A^{-1}$ $$ B\cdot [A,M_g]\cdot B = M_g \cdot B - B \cdot M_g = [M_g,B]=0 \label{8}\tag{8}$$
But I do also wish for precise bibliographical references and possible alternative arguments.
Such a statement is mentionned e.g. in "Symmetries and Laplacians, Introduction to Harmonic Analysis, Group Representations and Applications" by David Gurarie, §4.6 p.102 and it cites "Applications of Lie groups to differential equations" (GTM) by P. Olver, "Group analysis of differential equations" but L. V. Ovsiannikov and "Lie groups and special functions" and "Symmetry and separation of variables" by W. Miller, but I couldn't find what I expected and these are not so easy to read.