By definition, the marginal pdf of $X_1$ is given by :
$$ f_{X_1}(x) =\int_{x_2\in\mathbb R} f_{X_1,X_2}(x,x_2)\ d x_2 \quad\forall x\in\mathbb R\tag1$$
Where $f_{X_1,X_2}$ is the pdf of $(X_1,X_2)$. Of course, the marginal of $X_2$ is defined in an analogous way. Now if we write the joint distribution of $(X_1,X_2)$ as
$$\left(\begin{array}{c}X_1\\X_2 \end{array}\right) \sim \mathcal N \left[\left(\begin{array}{c}\mu_1\\ \mu_2 \end{array}\right), \left(\begin{array}{cc}\sigma^2_1 & \rho \sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma^2_2 \end{array}\right)\right] $$
we find that the joint pdf is given by
$$\begin{align}
&f_{X_1,X_2}(x_1,x_2)=\frac{1}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}} \times\\
&\exp \left\{-\frac{1}{2 (1-\rho^2)}\bigg[\bigg(\frac{x_1-\mu_1}{\sigma_1}\bigg)^2
+\bigg(\frac{x_2-\mu_2}{\sigma_2}\bigg)^2-2\rho \frac{(x_1-\mu_1)(x_2-\mu_2)}{\sigma_1 \sigma_2} \bigg] \right\}
\end{align}\tag2 $$
Furthermore, observe that the argument of the exponential above can be rewritten as
$$\begin{align}
&\bigg(\frac{x_1-\mu_1}{\sigma_1}\bigg)^2
+\bigg(\frac{x_2-\mu_2}{\sigma_2}\bigg)^2-2\rho \frac{(x_1-\mu_1)(x_2-\mu_2)}{\sigma_1 \sigma_2}\\
=&(1-\rho^2)\bigg(\frac{x_1-\mu_1}{\sigma_1}\bigg)^2 + \rho^2\bigg(\frac{x_1-\mu_1}{\sigma_1}\bigg)^2
+\bigg(\frac{x_2-\mu_2}{\sigma_2}\bigg)^2-2\rho \frac{(x_1-\mu_1)(x_2-\mu_2)}{\sigma_1 \sigma_2}\\
=& (1-\rho^2)\bigg(\frac{x_1-\mu_1}{\sigma_1}\bigg)^2 +\left(\frac{x_2 - \mu_2 - \rho\frac{\sigma_2}{\sigma_1}(x_1-\mu_1)}{\sigma_2}\right)^2 \end{align} $$
So plugging this back in $(2)$, we find that the joint pdf of $(X_1,X_2)$ is equal to
$$\begin{align}
&f_{X_1,X_2}(x_1,x_2)=\frac{1}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}} \times\\
&\exp \left[-\frac{1}{2} \hat{x_1}^2-\frac{1}{2}\left(\frac{x_2 - \mu_2 - \rho\sigma_2\hat{x_1}}{ \sigma_2\sqrt{1-\rho^2}}\right)^2 \right].
\end{align}\tag3 $$
Where we used the notation $\hat{x_1}:=(x_1-\mu_1)/\sigma_1$ for convenience. Finally, if we remember that
$$\int_{\mathbb R}e^{-(x-\mu)^2/2\sigma^2} dx= \sigma\sqrt{2\pi} \tag4$$
for all $\sigma,\mu$ and plug all our findings into $(1)$, we find that the marginal pdf of $X_1$ is given for all $x\in \mathbb R$ by
$$\begin{align}
f_{X_1}(x) &= \frac{e^{-\hat{x_1}^2/2}}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}} \cdot\int_{-\infty}^\infty \exp \left[-\frac{1}{2}\left(\frac{x_2 - \mu_2 - \rho\sigma_2\hat{x_1}}{ \sigma_2\sqrt{1-\rho^2}}\right)^2 \right] dx_2\\
&\stackrel{(4)}{=} \frac{e^{-\hat{x_1}^2/2}}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}}\cdot \sigma_2\sqrt{2\pi(1-\rho)^2}\\
&=\frac{e^{-\hat{x_1}^2/2}}{\sqrt{2 \pi} \sigma_1 }
\end{align}, $$
and finally, replacing $\hat{x_1}$, we find
$$ \boxed{f_{X_1}(x) = \frac{1}{\sqrt{2 \pi} \sigma_1}\cdot\exp\left(-\frac{1}{2}\left[\frac{x_1-\mu_1}{\sigma_1}\right]^2\right)}, $$
That is, the marginal distribution of $X_1$ is normal with mean $\mu_1$ and variance $\sigma_1^2$, regardless of the correlation between $X_1$ and $X_2$.
Note : if all you were interested in is proving that $X_1$ and $X_2$ are marginally Gaussian, this approach is needlessly tedious, and your original argument is the way to go. While the above method shows that the marginal distributions are (surprisingly ?) independent of the correlation $\rho$, this well known result applied to the linear transformations you mentioned could have done the same job for way less work !
