There's a nice discussion of local isometric immersions of hyperbolic planes in Robert McLachlan's expository article:
McLachlan, Robert, A gallery of constant-negative-curvature surfaces, Math. Intell. 16, No. 4, 31-37 (1994). ZBL0812.53003.
See also his accompanying website gallery and Richard Palais' gallery of pseudospherical surfaces. Ghys has a nice survey too in french, but I didn't see how to fill in the sketch of proof that he gives of Hilbert's theorem.
I'm not sure what Thurston had in mind referring to the twisting of the lines of curvature. Following Hilbert's proof, the asymptotic curves give natural coordinates on an isometrically immersed hyperbolic surface which form a Chebyshev net. In these coordinates, the metric (first fundamental form) has the form
$$\mathrm{I} = \left(\begin{array}{cc}1 & \cos(\omega) \\\cos(\omega) & 1\end{array}\right), $$ where $\omega$ is the angle between the asymptotic curves. Note that the area form is $dA=\sqrt{det(\mathrm{I})}dxdy=\sin(\omega)dxdy$.
A nice proof of this is given in the Proposition on p. 192 of Ghys' paper (he says, translated: "Let us try to demonstrate it by avoiding too many calculations. . ."). He shows that the shape operator acts on vectors in the two asymptotic directions by rotation by $\pi/2$ and $-\pi/2$. Then giving a local parameterization by coordinates parallel to the asymptotic directions, he uses this to show that the mixed second partial derivative of this parameterization is orthogonal to the surface, which is equivalent to the Chebyshev net condition. I'm not sure if there's an intuitive way to understand this. One thing to note is that the asymptotic curves have constant torsion ($\pm 1$), which I think is related to the result about the action of the shape operator.
McLachlan and Ghys describe a Chebyshev net as a piece of clothing with fibers going in two directions, where the fiber lengths and their crossings remain rigid, but the angles between them can vary. Then any such piece of clothing can be draped over a piece of a surface by varying the fiber angles, also like a fish net or fishnet stocking (indeed, Chebyshev considered the fitting of clothes in the presentation he made in 1878). See also Terng and Uhlenbeck "The Geometry of Solitons".
Then the second fundamental form has the formula $$\mathrm{I\!I} = \left(\begin{array}{cc}0 & \sin(\omega) \\\sin(\omega) & 0\end{array}\right). $$ The $0$ on the diagonal are equivalent to taking asymptotic curves for the coordinates. One checks that $\frac{det(\mathrm{I\!I})}{det(\mathrm{I})}+-1 = 0,$ agreeing with Gauss' Theorema Egregrium.
Now McLachlan notes that the Gauss-Codazzi-Mainardi equations imply that the angle $\omega$ satisfies the sine-Gordon equation and has $0 < \omega <\pi $, since the angle between the coordinates must be non-zero.
$$\frac{\partial^2\omega}{\partial x\partial y}= \sin(\omega).$$
On any rectangle in the coordinate chart, the sine-Gordon equation implies that the area is bounded by $2\pi$ ("Hazzidakis' Formula"). Let the rectangle $X=[x_1,x_2]\times[y_1,y_2]$, then
$$\omega(x_1,y_1)-\omega(x_1,y_2)-\omega(x_2,y_1)+\omega(x_2,y_2) = \iint_X \frac{\partial^2 \omega}{dxdy} dx dy \\ = \iint_X \sin(\omega) dx dy = \iint_{\psi(X)} dA = Area(X).$$
Since $0<\omega < \pi$, the area is bounded by $2\pi$.
Now suppose that a global isometric immersion exists $\varphi:\mathbb{H}^2\to \mathbb{E}^3$. Take a point $p$ in the hyperbolic plane, then we have canonical Chebyshev coordinates about $p$, since the asymptotic curves are complete (first follow one asymptotic line through $p$ a distance $x$, then follow the transverse asymptotic curve a distance $y$ to get to the coordinates $(x,y)$ to get $\psi:\mathbb{R}^2\to \mathbb{H}^2$). If we pull back the metric to this chart, then we cannot get a complete metric because the area is bounded by $2\pi$, but the area of the hyperbolic plane is infinite. Hence by Hopf-Rinow there must be some geodesic line $\gamma:\mathbb{R}\to \mathbb{R}^2$ in this chart that is finite length, so that $\underset{t\to ∞}{\lim} \psi(\gamma(t))=q\in \mathbb{H}^2$. But along this curve then we must have $\omega$ getting arbitrarily close to $0$ or $\pi$, since if $\omega$ remains bounded away from these, then the length of the tangent vector is bounded away from $0$, so the curve would have infinite length. The principal curvatures are $\cot(\omega)\pm \csc(\omega)=\frac{\cos(\omega)\pm 1}{\sin(\omega)}$, and thus as $\omega\to 0$ or $\omega \to \pi$, the principal curvatures get arbitrarily close to $0$ and $\infty$. But by completeness of the immersion, the limit point $\underset{t\to ∞}{\lim} \varphi(\psi(\gamma(t)))=\varphi(q)$ would have to exist with finite principal curvatures, a contradiction.