Let $f$ be an automorphism of the algebra of octonions. Is it true that $f$ preserves some quaternionic subalgebra? Has the statement an elementary proof?

Seen as a map of $8$-dimensional Euclidean vector spaces, $f$ is obviously (special) orthogonal, so we can find an orthonormal basis on which is has a block diagonal form of $2\times 2$ rotation matrices, and certainly at least $u,v$ unit octonions, orthogonal to each other and both orthogonal to $1$ (=pure imaginary) such that $f(u) = \cos\theta\cdot u + \sin\theta\cdot v$ and $f(v) = -\sin\theta\cdot u + \cos\theta\cdot v$. Then $1,u,v,uv$ span (and in fact, are an orthonormal basis of) a quaternion algebra which is preserved by $f$ (incidentally, the unit imaginary octonion $uv$ is also preserved by $f$).

I would start by finding imaginary octonion $i$ fixed by $f$. This we can conclude from knowing that $f$ is orthogonal map of $R^7$ perpendicular to $1$. Now we apply argument used in Gro-Tsen answer. Each element in $SO_7$ has fixed vector.

Next I would observe that multiplication by $i$ defines complex structure on perpendicular $R^6$. This complex structure is preserved by $f$, because $f(\color {red}ix)=f(\color{red}i)f(x)=\color{red}if(x)$.

We have proved that $f$ belongs to $U_3$. Now, there exists basis $\color{blue}u,\color{green}v,\color{brown}w$ in $\mathbb C^3$ such that $f$ is diagonal in it. This means $f(\color{blue}u)=e^{\color{red}i\alpha}\color{blue}u$. The subalgebra $\langle1,\color{red}i,\color{blue}u,\color{red}i\color{blue}u\rangle$ is quaternion subalgebra and it is fixed by $f$.

For example:

1) $\color{red}i\color{blue}u=-\color{blue}u\color{red}i$ because $\color{red}i,\color{blue}u$ are perpendicular and imaginary;

2) $(\color{red}i\color{blue}u)\color{blue}u=\color{red}i\color{blue}u^2=-\color{red}i$ (Moufang identity).

These are geometrical intuitions which I use in octonions. To have strict proof we need to have definition of octonion multiplication.

Side comment - interestingly for two perpendicular imaginary octonions $x,y$ element $(L_xR_y)^2$ is octonion automorphism. It is identity on quaternion subalgebra generated by $x, y$ and minus identity on perpendicular 4-space. Such elements form submanifold in $G_2$ being 8-dimensional exceptional symmetric space $G$ of all subalgebras $\{\mathbb H \subset \mathbb O\}$.

Second comment is that this formula $(L_xR_y)^2$ gives automorphism for octonions over finite fields. We need carefully define what "imaginary octonion" means in this case.