Those distributions you call "marginal" are not marginal. They are conditional distributions because you wrote $x \mid y$. The marginal distribution of $X$, for example, is necessarily independent of the value of $Y$.
To see how the conditional distribution is gamma, all you have to do is write $$f_{X \mid Y}(x) = \frac{f_{X,Y}(x,y)}{f_Y(y)} \propto f_{X,Y}(x,y).$$ That is to say, the conditional distribution is proportional to the joint distribution, appropriately normalized. So we have $$f_{X \mid Y}(x) \propto x^2 e^{-x(y^2+4)},$$ completely ignoring any factors that are not functions of $x$. Then we recognize that the gamma distribution has density $$f_S(s) \propto s^{a-1} e^{-bs},$$ so the choice of shape $a = 3$ and rate $b = y^2+4$ demonstrates that the conditional distribution $X \mid Y \sim \operatorname{Gamma}(a = 3, b = y^2+4)$.
The conditional distribution of $Y \mid X$ is done in a similar fashion. Just ignore constants of proportionality: $$f_{Y \mid X}(y) \propto e^{-(x+1)y^2+2y},$$ but this one requires us to complete the square to get it to look like a normal density: $$-(x+1)y^2 + 2y = (x+1)\left(-\left(y - \tfrac{1}{x+1}\right)^2 \right) + \tfrac{1}{x+1},$$ and after exponentiating and removing the $e^{1/(x+1)}$ factor, comparing this against $$f_W(w) \propto e^{-(w-\mu)^2/(2\sigma^2)},$$ we see that we have a normal density with mean $\mu = 1/(x+1)$ and variance $\sigma^2 = 1/(2(x+1))$.
Now, if you wanted the marginal distributions, you would need to integrate: $$f_X(x) = \int_{y=-\infty}^\infty f_{X,Y}(x,y) \, dy,$$ for example. And as you can see, this expression will not be a function of $Y$. The difference is that if I simulated realizations of ordered pairs $(X_i, Y_i)$ from the joint distribution, the marginal density for $X$ would be what you would see if I only told you the values of $X_i$. The conditional distribution of $X$ given $Y = y$ would be what you would see if I only gave you the $X_i$ for which the corresponding $Y_i$ was equal to $y$.
