There isn't any difference here between using curved spacetime and curvilinear coordinates. While one could tell "true" curvature in spacetime by, say, calculating scalars like $R^{\alpha\beta\gamma\delta} R_{\alpha\beta\gamma\delta}$ or $R^\mu{}_\mu$ and seeing if they vanish, the fact is your differential operator doesn't care. Put another way, all you care about is whether or not the connection coefficients $\Gamma^\sigma_{\mu\nu}$ vanish, not whether some special coordinate-independent quantity does.
Because of this, Minkowski spacetime often but not always implicitly involves Cartesian coordinates, otherwise one could probably just call it "flat." If your spacetime is flat but your coordinates are not Cartesian, there will be terms you omitted from your third equation. In any spacetime and for any scalar1 $\varphi$,
\begin{align}
\square\varphi & = g^{\mu\nu} \nabla_\mu \nabla_\nu \varphi \\
& = g^{\mu\nu} \nabla_\mu (\partial_\nu \varphi) \\
& = g^{\mu\nu} (\partial_\mu \partial_\nu \varphi - \Gamma^\sigma_{\mu\nu} \partial_\sigma \varphi) \\
& = (\partial_\mu \partial^\mu - g^{\mu\nu} \Gamma^\sigma_{\mu\nu} \partial_\sigma) \varphi. \tag{1}
\end{align}
So if in coordinates $(t,x,y,z)$ we have $g_{\mu\nu} = \mathrm{diag}(-1,1,1,1)$, then sure, $\Gamma^\sigma_{\mu\nu} = 0$ and $\square\varphi = \partial_\mu \partial^\mu \varphi$. The same flat spacetime in spherical coordinates $(t,r,\theta,\phi)$ -- $g_{\mu\nu} = \mathrm{diag}(-1,1,r^2,r^2\sin^2\theta)$ -- will however have some nonzero connection coefficients. In fact, you should recognize the second term in (1) as being necessary, given that there are single derivatives in the spherical Laplacian:
$$ \nabla^2\varphi = \left(\partial_r^2 + \frac{2}{r} \partial_r + \frac{1}{r^2} \partial_\theta^2 + \frac{1}{r^2\tan\theta} \partial_\theta + \frac{1}{r^2\sin^2\theta} \partial_\phi^2\right) \varphi. $$
The thing to remember is covariant derivatives reduce to partial derivatives generally only in the following two circumstances:
- $\nabla_\mu \varphi = \partial_\mu \varphi$ in any spacetime in any coordinates for a scalar $\varphi$.
- $\nabla_\mu T^\alpha{}_\beta = \partial_\mu T^\alpha{}_\beta$ in flat spacetime in Cartesian coordinates, where $g_{\mu\nu} = \eta_{\mu\nu}$, for any tensor $T$ (with as many indices as you want).
Other than those two times, you have to resort to the general rule
$$ \nabla_\mu T^\alpha{}_\beta = \partial_\mu T^\alpha{}_\beta + \Gamma^\alpha_{\mu\sigma} T^\sigma{}_\beta - \Gamma^\rho_{\mu\beta} T^\alpha{}_\rho, $$
where extra indices follow the same pattern as $\alpha$ or $\beta$ depending on if they are up or down.
1 Note that if $\varphi$ is a higher-rank tensor, the simplification $\nabla\varphi \to \partial\varphi$ no longer holds. Also, in the context of wave equations, often the differential operator you want is not the d'Alambertian but rather the de Rham operator, which adds a coupling between the Ricci tensor and your tensor being differentiated (a "curvature coupling"). This reduces to the d'Alambertian on scalars.
Gravity is a gauge theory. Gauge transformations are diffeomorphisms (coordinate changes) described by your equations. Therefore, the space of all possible metrics (the moduli space) is the quotient of the space of all $g_{\mu \nu}$ over these coordinate changes.
So your $g_{\mu \nu}$ can be set to $\eta_{\mu \nu}$ by some coordinate transformation. It means that they belong to the same equivalence class.
On the other hand, $q_{\mu \nu}$ belongs to another equivalence class. It can be seen by computing the Riemann curvature tensor. For any $g_{\mu \nu}$ it should be zero, but not for $q_{\mu \nu}$.
Best Answer
Mathematically speaking they are the same operator. Usually we reserve the d'Alembertian for 3+1 dimensional spacetime (so in absence of curvature it takes the form $\partial_0^2 - \nabla^2$), while the Laplace-Beltrami operator is defined for an aribtrary dimensional manifold with arbitrary signature. The only possible difference is that sometimes (not always, though), $\Box$ is defined as $\partial_0^2 - \partial_1^2 - \partial_2^2 - \partial_3^2$ independently of signature, so if your metric is $(-+++)$ then you will have $\Box = -\nabla^2$, where $\nabla^2$ here means the LB operator.