Another way to solve this is to change coordinates so the weighted inner product $\langle x,y\rangle_P$ turns into the ordinary dot product, and then do plain Gram-Schmidt in the new coordinates.
Since $P = \text{diag}(1,9,16)$, write $P = D^2$ where $D = \text{diag}(1,3,4)$. Then for any $x,y$, $\langle x,y\rangle_P = x^\top D^2 y = (Dx)\cdot(Dy)$, the ordinary dot product of the scaled vectors $u=Dx$ and $v=Dy$.
Scale the two spanning vectors using $D$: $u_1 = Dv_1 = \begin{pmatrix}1\\0\\0\end{pmatrix}$, and $u_2 = Dv_2 = \begin{pmatrix}1\\3\\4\end{pmatrix}$.
Now do ordinary Gram-Schmidt on $u_1, u_2$. The vector $u_1$ already has length 1, so $\hat e_1 = \begin{pmatrix}1\\0\\0\end{pmatrix}$.
The projection of $u_2$ on $\hat e_1$ is $(u_2\cdot \hat e_1)\hat e_1 = \begin{pmatrix}1\\0\\0\end{pmatrix}$, so the leftover part is $u_2 - \begin{pmatrix}1\\0\\0\end{pmatrix} = \begin{pmatrix}0\\3\\4\end{pmatrix}$, which has length $\sqrt{3^2+4^2}=5$.
So $\hat e_2 = \begin{pmatrix}0\\\frac{3}{5}\\\frac{4}{5}\end{pmatrix}$.
Convert back to the original coordinates using $x = D^{-1}u$, with $D^{-1}=\text{diag}(1,\frac{1}{3},\frac{1}{4})$. This gives $e_1 = D^{-1}\hat e_1 = \begin{pmatrix}1\\0\\0\end{pmatrix}$ and $e_2 = D^{-1}\hat e_2 = \begin{pmatrix}0\\\frac{1}{5}\\\frac{1}{5}\end{pmatrix}$, exactly matching the basis found by direct Gram-Schmidt with the $P$-inner product.
Let's summarize:
So the correct orthonormal basis is $\left\{\begin{pmatrix}1\\0\\0\end{pmatrix}, \begin{pmatrix}0\\\frac{1}{5}\\\frac{1}{5}\end{pmatrix}\right\}$, option (C).
\[ \boxed{\left\{\begin{pmatrix}1\\0\\0\end{pmatrix}, \begin{pmatrix}0\\\frac{1}{5}\\\frac{1}{5}\end{pmatrix}\right\}} \]