The $1/|r|^2$ dependence is fundamentally geometrical in nature, stemming from our 3d world.
Differential equations, even PDEs, are sometimes said to have Green's functions. Consider the differential equation
$$\nabla \cdot K = \alpha$$
for some vector field $K$ and some scalar field $\alpha$. This is structurally identical to Gauss's law. This equation has a Green's function $G$ that satisfies
$$\nabla \cdot G = \delta$$
where $\delta$ is the Dirac delta function. This describes, in essence, a point source $\delta$ generating a field $G$.
The Green's function is given by
$$G(r) = \frac{\hat r}{4\pi |r|^2}$$
Thus, any differential equation of the form $\nabla \cdot K = \alpha$ has the same Green's function in 3d--a point source always generates the same basic field.
In short, it is a physical statement to say $\nabla \cdot E = \rho/\epsilon_0$, but once that has been said, Coulomb's law, which describes a point charge, must inevitably follow due to the mathematical structure of the equations and the geometrical nature of 3d space.
(Of course, you're free to work the other way around.)
Edit: crucially, in 2d and 1d, the Green's functions are different. You probably already know these solutions. The 2d Green's function is something familiar from the field of a uniform line charge, and has $1/|r|$ dependence. The 1d Green's function has constant magnitude but changes direction on opposite sides of the "point" charge--this is the basic feature of a uniformly charged sheet.