I got stucked with the same problem and Zwähnia's solution is great and works perfectly in the case you works with Euler angle but couldn't help me with Tait-Bryan angle. There is kind of a bug and I couldn't track it down yet, so there is a workaround.
The problem I had comes from the fact you redefined the \tdplotcalctransformrotmain command (and so did I). This makes in some point tikz-3dplot go crazy.
As you can see in the first graph, using tikz-3dplot implementation of the rotation, the circle are drawn in the z_B - x_B plane (the blue circle), y_B - z_B (the red circle) and x_B - y_B (the black circle) whereas in the Tait-Bryan case the red and the black circle are in the same plane
The codes are similar though...
\documentclass{standalone}
\usepackage{tikz}
\usepackage{tikz-3dplot}
\usepackage{pgfplots}
% Workaround for making use of externalization possible
% -> remove hardcoded pdflatex and replace by lualatex
\usepgfplotslibrary{external}
% Redefine rotation sequence for tikz3d-plot to z-y-x
\newcommand{\tdseteulerxyz}{
\renewcommand{\tdplotcalctransformrotmain}{%
%perform some trig for the Euler transformation
\tdplotsinandcos{\sinalpha}{\cosalpha}{\tdplotalpha}
\tdplotsinandcos{\sinbeta}{\cosbeta}{\tdplotbeta}
\tdplotsinandcos{\singamma}{\cosgamma}{\tdplotgamma}
%
\tdplotmult{\sasb}{\sinalpha}{\sinbeta}
\tdplotmult{\sasg}{\sinalpha}{\singamma}
\tdplotmult{\sasbsg}{\sasb}{\singamma}
%
\tdplotmult{\sacb}{\sinalpha}{\cosbeta}
\tdplotmult{\sacg}{\sinalpha}{\cosgamma}
\tdplotmult{\sasbcg}{\sasb}{\cosgamma}
%
\tdplotmult{\casb}{\cosalpha}{\sinbeta}
\tdplotmult{\cacb}{\cosalpha}{\cosbeta}
\tdplotmult{\cacg}{\cosalpha}{\cosgamma}
\tdplotmult{\casg}{\cosalpha}{\singamma}
%
\tdplotmult{\cbsg}{\cosbeta}{\singamma}
\tdplotmult{\cbcg}{\cosbeta}{\cosgamma}
%
\tdplotmult{\casbsg}{\casb}{\singamma}
\tdplotmult{\casbcg}{\casb}{\cosgamma}
%
%determine rotation matrix elements for Euler transformation
\pgfmathsetmacro{\raaeul}{\cacb}
\pgfmathsetmacro{\rabeul}{\casbsg - \sacg}
\pgfmathsetmacro{\raceul}{\sasg + \casbcg}
\pgfmathsetmacro{\rbaeul}{\sacb}
\pgfmathsetmacro{\rbbeul}{\sasbsg + \cacg}
\pgfmathsetmacro{\rbceul}{\sasbcg - \casg}
\pgfmathsetmacro{\rcaeul}{-\sinbeta}
\pgfmathsetmacro{\rcbeul}{\cbsg}
\pgfmathsetmacro{\rcceul}{\cbcg}
}
}
% Set the plot display orientation
% Syntax: \tdplotsetdisplay{\theta_d}{\phi_d}
\tdplotsetmaincoords{60}{140}
\pgfmathsetmacro{\zRot}{10}
\pgfmathsetmacro{\yRot}{10}
\pgfmathsetmacro{\xRot}{10}
%%%%%%%%% Using standard euler angles implemented in default tikz-3dplot implementation
\begin{document}
%%%%%%%%%%%%% Z-Y-Z
\begin{tikzpicture}[scale=2.5,tdplot_main_coords]
% Set origin of main (body) coordinate system
\coordinate (O) at (0,0,0);
% Draw main coordinate system
\draw[red, ,->] (0,0,0) -- (1,0,0) node[anchor=north east]{$x_{\mathcal{I}}$};
\draw[red, ,->] (0,0,0) -- (0,1,0) node[anchor=north west]{$y_{\mathcal{I}}$};
\draw[red, ,->] (0,0,0) -- (0,0,1) node[anchor=south]{$z_{\mathcal{I}}$};
% Rotate to final frame
\tdplotsetrotatedcoords{\zRot}{\yRot}{\xRot}
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (1,0,0) node[anchor=west]{$x_{\mathcal{B}}$};
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (0,1,0) node[anchor=west]{$y_{\mathcal{B}}$};
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (0,0,1) node[anchor=south]{$z_{\mathcal{B}}$};
%Draws circle representing the rotated planes. Each of these should be "pointed" by two arrows.
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=black]{(0,0,0)}{1}{0}{350}{anchor=south west,color=black}{$x-y$}
\tdplotsetrotatedthetaplanecoords{0}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=blue]{(0,0,0)}{1}{0}{350}{anchor=south west,color=blue}{$x-z$}
\tdplotsetrotatedthetaplanecoords{90}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=red]{(0,0,0)}{1}{0}{350}{anchor=south west,color=red}{$y-z$}
\end{tikzpicture}
%%%%%% Change the rotation matrix in order to use Tait-Bryan angles
\tdseteulerxyz
%%%%%%%%%%%%% Z-Y-X
\begin{tikzpicture}[scale=2.5,tdplot_main_coords]
% Set origin of main (body) coordinate system
\coordinate (O) at (0,0,0);
% Draw main coordinate system
\draw[red, ,->] (0,0,0) -- (1,0,0) node[anchor=north east]{$x_{\mathcal{I}}$};
\draw[red, ,->] (0,0,0) -- (0,1,0) node[anchor=north west]{$y_{\mathcal{I}}$};
\draw[red, ,->] (0,0,0) -- (0,0,1) node[anchor=south]{$z_{\mathcal{I}}$};
% Rotate to final frame
\tdplotsetrotatedcoords{\zRot}{\yRot}{\xRot}
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (1,0,0) node[anchor=west]{$x_{\mathcal{B}}$};
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (0,1,0) node[anchor=west]{$y_{\mathcal{B}}$};
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (0,0,1) node[anchor=south]{$z_{\mathcal{B}}$};
%Draws circle representing the rotated planes. Each of these should be "pointed" by two arrows.
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=black]{(0,0,0)}{1}{0}{350}{anchor=south west,color=black}{$x-y$}
\tdplotsetrotatedthetaplanecoords{0}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=blue]{(0,0,0)}{1}{0}{350}{anchor=south west,color=blue}{$x-z$}
\tdplotsetrotatedthetaplanecoords{90}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=red]{(0,0,0)}{1}{0}{350}{anchor=south west,color=red}{$y-z$}
\end{tikzpicture}
\end{document}
I wanted to preserve the use of the angle macros \xRot, \yRot and \zRot which is not possible anymore if you set up the tdplotsetrotatedthetaplanecoords as does Zwähnia.
My solution consists in turning the coords backwards to find easily the good theta plane.
\documentclass{standalone}
\usepackage{tikz}
\usepackage{tikz-3dplot}
\usepackage{pgfplots}
% Workaround for making use of externalization possible
% -> remove hardcoded pdflatex and replace by lualatex
\usepgfplotslibrary{external}
% Redefine rotation sequence for tikz3d-plot to z-y-x
\newcommand{\tdseteulerxyz}{
\renewcommand{\tdplotcalctransformrotmain}{%
%perform some trig for the Euler transformation
\tdplotsinandcos{\sinalpha}{\cosalpha}{\tdplotalpha}
\tdplotsinandcos{\sinbeta}{\cosbeta}{\tdplotbeta}
\tdplotsinandcos{\singamma}{\cosgamma}{\tdplotgamma}
%
\tdplotmult{\sasb}{\sinalpha}{\sinbeta}
\tdplotmult{\sasg}{\sinalpha}{\singamma}
\tdplotmult{\sasbsg}{\sasb}{\singamma}
%
\tdplotmult{\sacb}{\sinalpha}{\cosbeta}
\tdplotmult{\sacg}{\sinalpha}{\cosgamma}
\tdplotmult{\sasbcg}{\sasb}{\cosgamma}
%
\tdplotmult{\casb}{\cosalpha}{\sinbeta}
\tdplotmult{\cacb}{\cosalpha}{\cosbeta}
\tdplotmult{\cacg}{\cosalpha}{\cosgamma}
\tdplotmult{\casg}{\cosalpha}{\singamma}
%
\tdplotmult{\cbsg}{\cosbeta}{\singamma}
\tdplotmult{\cbcg}{\cosbeta}{\cosgamma}
%
\tdplotmult{\casbsg}{\casb}{\singamma}
\tdplotmult{\casbcg}{\casb}{\cosgamma}
%
%determine rotation matrix elements for Euler transformation
\pgfmathsetmacro{\raaeul}{\cacb}
\pgfmathsetmacro{\rabeul}{\casbsg - \sacg}
\pgfmathsetmacro{\raceul}{\sasg + \casbcg}
\pgfmathsetmacro{\rbaeul}{\sacb}
\pgfmathsetmacro{\rbbeul}{\sasbsg + \cacg}
\pgfmathsetmacro{\rbceul}{\sasbcg - \casg}
\pgfmathsetmacro{\rcaeul}{-\sinbeta}
\pgfmathsetmacro{\rcbeul}{\cbsg}
\pgfmathsetmacro{\rcceul}{\cbcg}
}
}
% Set the plot display orientation
% Syntax: \tdplotsetdisplay{\theta_d}{\phi_d}
\tdplotsetmaincoords{60}{140}
\pgfmathsetmacro{\zRot}{10}
\pgfmathsetmacro{\yRot}{10}
\pgfmathsetmacro{\xRot}{10}
%%%%%%%%% Using standard euler angles implemented in default tikz-3dplot implementation
\begin{document}
%%%%%% Change the rotation matrix in order to use Tait-Bryan angles
\tdseteulerxyz
%%%%%%%%%%%%% Z-Y-X
\begin{tikzpicture}[scale=2.5,tdplot_main_coords]
% Set origin of main (body) coordinate system
\coordinate (O) at (0,0,0);
% Draw main coordinate system
\draw[red, ,->] (0,0,0) -- (1,0,0) node[anchor=north east]{$x_{\mathcal{I}}$};
\draw[red, ,->] (0,0,0) -- (0,1,0) node[anchor=north west]{$y_{\mathcal{I}}$};
\draw[red, ,->] (0,0,0) -- (0,0,1) node[anchor=south]{$z_{\mathcal{I}}$};
% Rotate to final frame
\tdplotsetrotatedcoords{\zRot}{\yRot}{\xRot}
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (1,0,0) node[anchor=west]{$x_{\mathcal{B}}$};
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (0,1,0) node[anchor=west]{$y_{\mathcal{B}}$};
\draw[thick,tdplot_rotated_coords,->, cyan] (0,0,0) -- (0,0,1) node[anchor=south]{$z_{\mathcal{B}}$};
%Draws circle representing the rotated planes. Each of these should be "pointed" by two arrows.
%\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=black]{(0,0,0)}{1}{0}{350}{anchor=south west,color=black}{$x-y$}
\tdplotsetrotatedthetaplanecoords{0}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=blue]{(0,0,0)}{1}{0}{350}{anchor=south west,color=blue}{$x-z$}
\tdplotsetrotatedthetaplanecoords{90}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=black]{(0,0,0)}{1}{0}{350}{anchor=south west,color=black}{$y-z$}
% Instead, rotates in the opposite direction
\tdplotsetrotatedcoords{\zRot+90}{\yRot}{\xRot}
%Draw the last circle in the right plane
\tdplotsetrotatedthetaplanecoords{0}
\tdplotdrawarc[dashed,tdplot_rotated_coords,->,color=red]{(0,0,0)}{1}{0}{350}{anchor=south west,color=red}{$x-z$}
\end{tikzpicture}
\end{document}
It's not really pretty but it has the advantage of staying with the use of the angle macros and allows to draw the sequence as I wanted it (I finish my last tikzpicture and I update my answer with it).
Hope this helps.
Cheers.
\tkzMarkAnglefromtkz-euclidewould be the right way to go (see TikZ: Create right angle triangle with angle labels), but in don't know if with some small adjstment could work also in 3D case.