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?

Here is another counterexample which gives closed-form distributions of $X, Y, X|Y = y$ and $Y|X = x$.
Let $Y, Z \text{ i.i.d. } \sim N(0, 1)$, and define $X = \frac{Z}{Y}$.  Then for $y \neq 0$ (the probability of $Y = 0$ is zero),
\begin{align}
X | Y = y \sim N(0, y^{-2}).
\end{align}
On the other hand, it is well known that the marginal distribution of $X$ is Cauchy distribution, i.e.,
\begin{align}
f_X(x) = \frac{1}{\pi(1 + x^2)}, \quad x \in \mathbb{R}. \tag{1}
\end{align}
And the joint distribution of $X$ and $Y$ can be evaluated as (where $\Phi$ and $\varphi$ denote CDF and PDF of the standard normal distribution respectively):
\begin{align}
 & F(x, y) = P[X \leq x, Y \leq y] \\
=& P[Z \leq Yx, Y \leq y, Y > 0] + P[Z \geq Yx, Y \leq y, Y < 0] \\
=& \begin{cases}
\int_{-\infty}^y(1 - \Phi(tx))\varphi(t)dt & y < 0, \\[1em]
\int_0^y\Phi(tx)\varphi(t)dt & y > 0. 
\end{cases}
\end{align}
Therefore, the joint density of $(X, Y)$ is given by
\begin{align}
 & f(x, y) = \frac{\partial^2F(x, y)}{\partial x\partial y} \\
=&\begin{cases}
-y\varphi(y)\varphi(yx) & y < 0, \\[1em]
y\varphi(y)\varphi(yx) & y > 0
\end{cases} \\
=& \frac{1}{2\pi}|y|e^{-(1 + x^2)y^2/2}. 
\tag{2}
\end{align}
$(1)$ and $(2)$ together yield the conditional density of $Y$ given $X = x$:
\begin{align}
f_{Y|X}(y|X = x) = \frac{f(x, y)}{f_X(x)} = 
\frac{1}{2}|y|(1 + x^2)e^{-(1 + x^2)y^2/2}. \tag{3}
\end{align}
Obviously, $(3)$ is not the density of any normal distribution (with $y$ as the variate).  Thus $Y | X = x$ is not normal.  For example, when $x = 0$, $(3)$ looks like as follows:

P.S., PDF $(3)$ may be termed as "double generalized gamma distribution", based on these two articles: Generalized gamma distribution and Double Gamma Distribution. The parameters linked to the generalized gamma distribution are $a = \sqrt{2(1 + x^2)^{-1}}$ (scale) and $d = 2, p = 2$.
