I have opened two textbooks and found the same notation in both. The notation for the center of mass in a lamina given by a variable density function $\rho(x,y)$ is $(\bar{x}, \bar{y}) = \left( \dfrac{M_y}{M}, \dfrac{M_x}{M}\right)$. This is where $$M =\iint_R \rho(x,y) dA$$ $$M_x =\iint_R y\rho(x,y) dA$$ $$M_y =\iint_R x\rho(x,y) dA$$
Why is the notation $M_x$ used when we are referring to the y-axis and vice versa?
The moment about the $y$-axis depends on the distance from the $y$-axis. The distance of a point $(x,y)$ from the $y$-axis is $x$. The integral represented by $M_x$ has the $x$ because one integrates the function against $x$, which is, again, the distance from the $y$-axis.