Suppose $X \mid Y$ and $Y$ are normally distributed. Does it follow that $Y \mid X$ is normally distributed? Suppose $X \mid Y$ and $Y$ are normally distributed. Does it follow that $Y \mid X$ is normally distributed?
And if so, how would one prove this?

whuber shows, by means of a counterexample, that the product of a Gaussian $X|Y$ times a Gaussian r.v. $Y$ does not necessarily lead to a joint Gaussian distribution, which in this case certainly doesn't have a conditional Gaussian.
Here is a (numerical) counter-example for your claim. Let $Y|X=x \sim N(x^3, 1)$ and $X \sim N(0,1)$. Then $(X,Y)$ is not jointly normal. The marginal of $Y$ is
$$
f_Y(y) = \frac{1}{2\pi}\int_{-\infty}^{\infty} e^{-\frac{1}{2}(y - x^3)^2 - \frac{1}{2}x^2}dx.
$$
Unfortunately this integral (as far as I know) cannot be computed analytically, but we can provide a fairly good approximation via adaptive quadratures. Thus the conditional we require is
$$
f_{X|Y} = \frac{e^{-\frac{1}{2}(y - x^3)^2 - \frac{1}{2}x^2}}{\int_{-\infty}^{\infty} e^{-\frac{1}{2}(y - x^3)^2 - \frac{1}{2}x^2}dx}.
$$
The result is disproved provided we are able to show that this conditional is not Gaussian. Toward this aim, let's fix $y = -1$. An elegant proof may be given by studying the properties of this function but here  I am taking a brute-force approach in which I approximate this beast numerically. If the conditional is gaussian we expect it to be unimodal and symmetric. This is the contradiction I'll be looking for.
 # propto the joint distribution
f_yx <- function(y, x) {
  exp(-0.5*(y-x^3)^2 - 0.5*x^2)
}

# the marginal of Y
f_y <- function(y) {
  integrate(function(t) f_yx(y, t), lower = -Inf, Inf)$value
}

# the conditional of x given y
x_given_y = function(x, y) f_yx(y,x)/f_y(y) 

# fix y
y = -1


x <- seq(-3, 3, len=100)
cond_y <- sapply(x, x_given_y, y=y)
plot(x, cond_y, type ="l", lwd=2)


If we can trust numerical integration (I'm using the integrate function here which is extremely robust!), we can see from the plot that this conditional density is definitely non-Gaussian. This contradicts the claim.
Side comment. There exist non-Gaussian bivariate distributions which have Gaussian conditional densities. One example of this is
$$
f(x,y) = C\exp\left(-(1+x^2)(1+y^2)\right),\quad -\infty <x,y<\infty,$$
where $C$ is the normalising constant. You can check that both $Y|X=x$ and $X|Y=x$ are Gaussian with suitable parameters.
