It is a good idea (not too natural, but not too out of nowhere) to consider $z = x + y + c + 1$ for some $c\in [2y,2x]$. (The $x+y$ choice is because of symmetry between the variables and is of the same order as the endpoints of the interval in the problem condition and the $+1$ could be thought just to relax the inequality condition for $c$ or could be thought much later after some manipulations.) Suppose, for contradiction, that $d>1$ divides $x+y+z = 2x+2y+c+1$ and $x^2 + xy + y^2$. Doing the clever thing of "conjugating the $c$", we get that $d$ divides $(2x+2y+1+c)(2x+2y+1-c) = 4(x^2 + y^2 + xy) + 4(xy + x + y) - c^2 + 1$. But with $xy+x+y = a^2$ and $c = 2a$ we get that $d$ must divide $1$, contradiction! (It is easy to check that indeed $2\sqrt{xy+x+y} \in [2y, 2x]$.)

The main idea is to find a positive integer $z$ such that $xy+yz+zx$ and $x+y+z+xy+yz+zx$ are both perfect square. This is related to an old problem from Kvant journal, i.e., Kvant 1799, stating if $xy+y+x$ is a perfect square then there is a positive integer $z$ such that $$x+yz, y+xz, z+xy, yz+y+z, zx+z+x, xy+yz+zx, x+y+z+xy+yz+zx,$$are all perfect square. Indeed, taking $t=\sqrt{x+y+xy}$ and $z=x+y+2t+1$ works. It also works in the problem at hand.