By using $x + y + z = 1$, we have
$$1+xy+yz+zx = (x+y+z)^2+xy+yz+zx = (x+y)(y+z)+(y+z)(z+x)+(z+x)(x+y).$$
With this equation and Cauchy-Schwarz inequality, we can deduce that
$$\begin{aligned} (LHS) &= \frac{1}{9} \left( \frac{1}{1-x} + \frac{1}{1-y} + \frac{1}{1-z} \right) \left( (x+3x^3) + (y+3y^3) + (z+3z^3) \right) \\ &\ge \left( \sqrt{\frac{3x^3+x}{9(1-x)}} + \sqrt{\frac{3y^3+y}{9(1-y)}} + \sqrt{\frac{3z^3+z}{9(1-z)}} \right)^2 \end{aligned}$$
Therefore, it is enough to show that for any real number $s \in (0, 1)$, the inequality
$$\frac{3s^3 + s}{9(1-s)} \ge \left( \frac{s\sqrt{1+s}}{\sqrt[4]{3+9s^2}} \right)^2$$
holds. If we expand the above inequality, then it is easy to check that this is equivalent to $3(9s^2 - 1)^2 \ge 0$, thus solving the problem. From the last inequality, we can check that the equality holds when $x = y = z = \frac{1}{3}$.
$\square$