We start by characterizing monomials of $R(x^2y, y^2z, z^2x)$. Claim: Polynomial $T(x,y,z)$ can be written as $R(x^2y, y^2z, z^2x)$ if for each its monomial $x^i y^j z^k$, $9|4i-2j+k$. Proof: We essentially need existance of $a,b,c \in \mathbb{Z}^{+}$ so that $i = 2a + b$, $j = 2b + c$, $k = 2c + a$. We can set $a = \frac{4i-2j+k}{9}$, $b = \frac{4j-2k+i}{9}$, $c = \frac{4k-2i+j}{9}$. We can check that one of the divisibilities implies others. We also need that all of them are non-negative, which we will handle later. Now we will be done if $P(x,y,z)Q(x,y,z) = T(x^9,y^9,z^9)$ then the condition of the last claim holds. For this, we use the following variation of the problem. Claim: Prove that for any non-zero polynomial $P\in \mathcal{M}$ there exists non-zero polynomials $Q,R\in \mathcal{M}$ such that\[ R(x^3,y,z) = P(x,y,z)Q(x,y,z). \] Proof: For a polynomial $P(x,y,z)$ and $i = 0,1,2$ , let $P_i(x,y,z)$ be the sum of all its monomials which have their $x$ power congruent to $i \pmod{3}$. Then $P(x,y,z) = P_0(x,y,z) + P_1(x,y,z) + P_2(x,y,z)$ and set $Q(x,y,z) = P_0(x,y,z)^2 + P_1(x,y,z)^2 + P_2(x,y,z)^2 - P_0(x,y,z)P_1(x,y,z) - P_1(x,y,z)P_2(x,y,z) - P_2(x,y,z)P_0(x,y,z)$ then $P(x,y,z)Q(x,y,z) = P_0(x,y,z)^3 + P_1(x,y,z)^3 + P_2(x,y,z)^3 - 3P_0(x,y,z)P_1(x,y,z)P_2(x,y,z)$. Each monomial is composed either of a product of three monomials from the same $P_i$ or three from different, both of which have $x$ power divisible by three. So every monomial of $PQ$ has $x$ power divisible by three and we can say that $P(x,y,z)Q(x,y,z) = R(x^3,y,z)$. Alternatively, one can check that $P(x,y,z)P(\omega x, y,z)P(\omega^2x, y , z) = R(x^3, y, z)$, where $\omega^2 + \omega + 1 =0$ and $R \in \mathcal{M}$, and so $P(\omega x, y,z)P(\omega^2x, y , z) = Q(x,y,z) \in \mathcal{M}$ too and it works. Now just spam the claim, interchanging variables, to get $R_1(x^3,y,z)$, then $R_2(x^9, y, z)$, then $R_3(x^9,y^3,z)$, $\dots$, $T(x^9,y^9,z^9)$. Every time we multiplied by polynomials from $\mathcal{M}$; so their product, namely $\frac{T(x^9,y^9,z^9)}{P(x,y,z)}$, is also in that set, and it works. Now, we can multiply both sides by $(xyz)^{3N}$ for a big enough $N$ to ensure that all the degrees are non-negative.

Using the notation from last solution, let $P_i(x,y,z)$ for $i = 0,1, \ldots, 8$ be the sum of monomials $x^a y^b z^c$ of $P$ that have $4a-2b+c \equiv i \pmod{9}$. Let's say that the monomials $x^a y^b z^c$ for which $4a-2b+c \equiv 0 \pmod{9}$ form a subset $\mathbb{S}$ and say that $\mathbb{S}[x,y,z]$ are the polynomials, all of which monomials are in $\mathbb{S}$ and have rational coefficients. For convinience, $\mathbb{S}[x,y,z] = \mathcal{S}$. We want to choose $Q(x,y,z) \in \mathcal{M}$ so that $P(x,y,z)Q(x,y,z) \in \mathcal{S}$. Suppose that $Q(x,y,z) = Q_0(x,y,z) + Q_1(x,y,z) + \dots + Q_8(x,y,z)$. Then what condition asks us is \[ \forall 1 \leq r \leq 8, \sum_{i+j \equiv r \pmod{9}}^{} P_i(x,y,z)Q_j(x,y,z) = 0 \]Which is the same as saying that \[ \forall 1 \leq r \leq 8, \sum_{i+j \equiv r \pmod{9}}^{} P_i(x,y,z)z^kQ_j(x,y,z)z^s= 0 \]where $k = -i \pmod{9}, s = -j \pmod{9}$, $0 \leq s \leq 8$ and $0<k+s<9$. Then $k+s$ is fixed and the equation holds. Note that k can be negative. The advantage of this is that now we can substitute $Q_i(x,y,z)z^s=Q'_i(x,y,z)$ where each and $Q'$ are in $\mathcal{S}$. For convenience, multiply each equation by $z^{9m}$ for large enough $m$ so that $P_i(x,y,z)z^{k+9m}$ is in $\mathcal{S}$ too. Now we have a homogeneous system of 8 linear equations with 9 variables, in which coefficients and variables are polynomials in $\mathcal{S}$. Now, as with usual systems of equations, we can get a non-trivial solution that can be represented as a rational function in polynomials from $\mathcal{S}$. But we can multiply the solutions by any polynomial. So just multiply it so that every $Q'$ is, in fact, polynomial. One remarkable property of $\mathcal{S}$ is that whenever we multiply two elements of it, we get an element of it, same with summation. So after we get a solution in polynomials for $Q'$, we know that they all are from $\mathcal{S}$. To get a solution in $Q$'s, we just need to divide by some $z$ in powers less than 9, we can ensure that we get polynomials by multiplying the solution for $Q'$ by a large power of $z$ first and then performing the division. This way we get the desired $Q(x,y,z)$. As in the previous, we are left to multiply both sides by $(xyz)^{3N}$ for a big enough $N$ to sort out negative powers.

Just let $x^iy^jz^k = (x^2y)^{a}(y^2z)^{b}(z^2x)^{c}$ and the iff condition follows, for the other direction we get a system of equations which clearly have solutions. (2022 C4 anyone?)

I had first gone with the most obvious idea of symmetric polynomials, however to show that if $P,Q$ can be written in that way, then $P+Q$ can also be written went hard, next the only way left to proceed is to find a characterization of monomials that can be represented in the said form, which leads to a solution!