If we express the change in anglular velocity as $\Delta\vec\omega$ in local coordinates, with for example the angles $\phi$, $\theta$ and $\psi$ rotations about the $Z$, $X$ and $Z$ respectively then the answer is
$$\Delta\vec{\omega}=\dot{\phi}\hat{k}+\mathrm{{Rot}(\hat{k},\phi})\left(\dot{\theta}\hat{i}+\mathrm{{Rot}(\hat{i},\theta})\left(\dot{\psi}\hat{k}\right)\right)$$
Where Rot(axis,angle) is a 3x3 rotation matrix. We apply the angular rotation components in sequence on the local coordinate axes this way. Also $\hat{i}=(1,0,0)$, $\hat{j}=(0,1,0)$ and $\hat{k}=(0,0,1)$.
Line of node must be the common normal between the two $z$-axes. We typically denote that as the $x$-axis (see Denavit-Hartenberg notation).
Cool question.
Let's try to formalize the general concepts presented in the other two responses thus far. Intuitively, the main idea is that the instantaneous axis of rotation is determined by the rotation that takes the object from its configuration at one instant to its configuration at the next instant as opposed to the rotation that takes it from its initial configuration to its configuration at some later time. Let's try to make this precise.
As you already wrote, for any time $t$ we have the following for the position $\vec x(t)$ of a point in the body at time $t$:
\begin{align}
\vec x(t) = R(t)\vec x(0) \tag{1}
\end{align}
where $R(t)$ is a rotation for each $t$. Now consider the configuration of the object at a neighboring time $t+\epsilon$. How is the position of a point in the body at this later point related to its position at the point $t$? Well, notice that
\begin{align}
\vec x(t+\epsilon)
&= R(t+\epsilon)\vec x(0) \\
&= \Big(R(t) + \dot R(t)\epsilon\Big)\vec x(0) + O(\epsilon^2) \\
&= \Big(I + \epsilon\dot R(t)R(t)^t\Big)\vec x(t) + O(\epsilon^2)
\end{align}
For small $\epsilon$, the order $\epsilon^2$ terms become insignificant, so the transformation that takes $\vec x(t)$ to $\vec x(t+\epsilon)$ is the gadget in big parentheses. Let's give it a name:
\begin{align}
\rho(t, \epsilon) = I+\epsilon\dot R(t)R(r)^t
\end{align}
But what exactly is this guy? Well, our intuition tells us that it should be a rotation with axis $\vec\omega(t)$; in particular, $\vec\omega(t)$ should be an eigenvector of this transformation. To see that this is the case,
I make the following:
Claim. $\rho(t,\epsilon)$ takes precisely the form of a "small" rotation about the axis $\vec\omega(t)$. In particular, $\vec \omega(t)$ is an eigenvector of $\rho(t,\epsilon)$.
Proof. First, recall that omega can be defined as the vector that generates the motion of the points in the body as follows:
\begin{align}
\dot{\vec x}(t) = \vec\omega(t)\times\vec x(t) \tag{2}
\end{align}
But notice, using eq. $(1)$ that
\begin{align}
\dot{\vec x}(t) = \dot R(t) \vec x(0) = \dot R(t)R(t)^t\vec x(t) \tag{3}
\end{align}
Now we notice that the matrix $\dot R(t) R(t)^t$ is antisymmetric since $R(t)$ is an orthogonal matrix for each $t$ (I'll leave this for you to prove. It follows that there exist functions $w_x(t), w_y(t), w_z(t)$ for which
\begin{align}
\dot R(t)R(t)^t = \begin{pmatrix}
0 & -w_z(t) & w_y(t) \\
w_z(t) & 0 & -w_x(t) \\
-w_y(t) & w_x(t) & 0 \\
\end{pmatrix}
\end{align}
Using this expression and comparing $(2)$ and $(3)$, it is not hard to show that $w_x,w_y,x_z$ are precisely the components of the angular velocity $\vec\omega$. Putting this all together, we see that $\rho(t,\epsilon)$ can be written as follows:
\begin{align}
\rho(t,\epsilon) = I + \epsilon\vec\omega(t)\times
\end{align}
where $\vec\omega(t)\times$ is shorthand for the transformation mapping any vector $\vec v$ to $\vec\omega(t)\times v$. Now, it is easy to show that $\vec\omega(t)$ is an eigenvector of $\rho(t,\epsilon)$;
\begin{align}
\rho(t,\epsilon)\vec\omega(t) = (I+\epsilon \vec\omega(t)\times)\omega(t) = \vec\omega(t) + \epsilon\vec\omega(t)\times\vec\omega(t) = \vec\omega(t).
\end{align}
It remains to show that $\rho(t,\epsilon)$ has the form of a "small" rotation about the axis $\vec\omega(t)$. Well, what does such a rotation look like? Well, a general rotation by an angle $\theta$ about an axis defined by a unit vector $\hat n$ can be written as follows (I'll again leave it to you to either show this, or find out why this is true):
\begin{align}
R(\hat n, \theta) = \exp(\theta\hat n\cdot \vec J)
\end{align}
Where $\vec J = (J_x, J_y, J_z)$ is a vector of matrices, the so-called rotation generators, whose components are defined by
\begin{align}
(J_i)_{jk}= -\epsilon_{ijk}
\end{align}
and $\exp$ is the matrix exponential. Notice that to first order in $\theta$, this shows that
\begin{align}
R(\hat n, \theta) = I + \theta\hat n\cdot J_i
\end{align}
But now notice that if we apply this to any vector $\vec v$, then we obtain
\begin{align}
R(\hat n, \theta)\vec v
&= (I+\theta\hat n\cdot \vec J)\vec v \\
&= \vec v + \theta n_i (J_i)_{jk}\vec v_k \hat e_j \\
&= \vec v - \theta n_i\epsilon_{ijk} v_k \hat e_j \\
&= \vec v + \theta \hat n\times \vec v \\
&= (I + \theta\hat n\times )\vec v
\end{align}
We now see by inspection that $\rho(t,\epsilon)$ has precisely the form of the first order approximation to a rotation $R(\hat n, \theta)$ if we simply take $\hat n = \hat \omega(t)$ and $\theta = \epsilon|\vec\omega(t)|$. $\blacksquare$
Best Answer
The theorem does not say that the actual axis of rotation is fixed. It says that the final configuration can be obtained by a rotation about a single axis. For instance, think about a sphere with its center fixed. Imagine the most general finite motion of this sphere. It will be a composition of many small rotations about different axis. Equivalently, the final configuration can be obtained by a single finite rotation about a fixed axis through its center. This was expected if we recall that rotations form a group. Thus the composition of several rotations is equal to a single rotation.