We shall prove the identity \begin{equation} \sum_{k=0}^n16^k[6\tan^4(2^kx)+8\tan^2(2^kx)+2]=16^{n+1}[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]- [6\cot^4(x)+8\cot^2(x)+2] \end{equation} via complex analysis.

Let $\mathcal{S}=3S_1+2S_2+2S_3$ denote the LHS of the above equality, and notice that $\mathcal{S}$ is meromorphic, and $\pi$-periodic. Note that the poles of $\tan^4(x)$ and $\tan^2(x)$ lie at points of the form $x_0=(m+\frac12)\pi$, with corresponding Laurent expansions \begin{equation} \begin{split} \tan^4(x)&=\frac1{(x-x_0)^4}-\frac4{3(x-x_0)^2}+O(1)\\ \tan^2(x)&=\frac1{(x-x_0)^2}+O(1)\\ \end{split} \end{equation} about $x_0$. Therefore the set of poles of $\mathcal{S}$ is given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole may be written uniquely as $x_0=\frac{(2m+1)\pi}{2^{k+1}}$ for integer $m$, and $0\leq k\leq n$, with corresponding Laurent expansion \begin{equation} \begin{split} \mathcal{S}(x)&=6\cdot 16^k\left[\frac1{2^{4k}(x-x_0)^4}-\frac4{3\cdot 2^{2k}(x-x_0)^2}\right] -8\cdot 16^k\left[\frac1{2^{2k}(x-x_0)^2}\right]+O(1)\\ &=\frac{6}{(x-x_0)^4}+O(1) \end{split} \end{equation} about $x_0$. Using the same logic as above, we may construct the meromorphic, $\pi$-periodic function $$f(x)=16^{n+1}[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]-[6\cot^4(x)+8\cot^2(x)+2]$$ which also has a set of poles given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole $x_0$ also has the corresponding Laurent expansion $f(x)=\frac{6}{(x-x_0)^4}+O(1)$ about $x_0$.

Finally, notice that $\lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\tan^2(x)=-1$, so \begin{equation} \lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\mathcal{S}(x)-f(x)=\sum_{k=0}^n16^k[6-8+2]-16^{n+1}[6-8+2]+[6-8+2]=0 \end{equation}

Since $\mathcal{S}$ and $f$ are both meromorphic and $\pi$-periodic, and share the same poles with the same Laurent expansions about those poles up to $O(1)$, then $\mathcal{S}-f$ is analytic and $\pi$-periodic, and since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $S-f$ is bounded. Therefore by Louville's theorem, $\mathcal{S}-f$ is constant. Finally, since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $\mathcal{S}-f$ is precisely $0$, which completes the proof.

Addendum: Note that this procedure generalises, and can be used to construct and prove higher order identities of this type. In fact, we have the following theorem which may be proven in the same way as above:

Let $c_{k,\ell}$ denote the $(-2k)$-th coefficient of the Laurent expansion of $\cot^{2\ell}(x)$ about $0$. If $a_1,\ldots, a_m$ is the solution to \begin{equation} \left[\begin{matrix} 0\\0\\0\\\vdots\\0\\1\\ \end{matrix}\right] = \left[\begin{matrix} 1&-1&1&\cdots&(-1)^{m-1}&(-1)^m\\ 0&c_{1,1}&c_{1,2}&\cdots&c_{1,m-1}&c_{1,m}\\ 0&0&c_{2,2}&\cdots&c_{2,m-1}&c_{2,m}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&c_{m-1,m-1}&c_{m-1,m}\\ 0&0&0&\cdots&0&c_{m,m}\\ \end{matrix}\right] \left[\begin{matrix} a_0\\a_1\\a_2\\\vdots\\a_{m-1}\\a_m\\ \end{matrix}\right] \end{equation} then \begin{equation} \sum_{k=0}^n2^{2mk}\sum_{\ell=0}^ma_\ell\tan^{2\ell}(2^k x)=\sum_{\ell=0}^m a_\ell[2^{2m(n+1)}\cot^{2\ell}(2^{n-1}x)-\cot^{2\ell}(x)] \end{equation}

We shall prove the identity \begin{equation} \sum_{k=0}^n16^k[6\tan^4(2^kx)+8\tan^2(2^kx)+2]=16^{n+1}[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]- [6\cot^4(x)+8\cot^2(x)+2] \end{equation} via complex analysis.

Let $\mathcal{S}=3S_1+2S_2+2S_3$ denote the LHS of the above equality, and notice that $\mathcal{S}$ is meromorphic, and $\pi$-periodic. Note that the poles of $\tan^4(x)$ and $\tan^2(x)$ lie at points of the form $x_0=(m+\frac12)\pi$, with corresponding Laurent expansions \begin{equation} \begin{split} \tan^4(x)&=\frac1{(x-x_0)^4}-\frac4{3(x-x_0)^2}+O(1)\\ \tan^2(x)&=\frac1{(x-x_0)^2}+O(1)\\ \end{split} \end{equation} about $x_0$. Therefore the set of poles of $\mathcal{S}$ is given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole may be written uniquely as $x_0=\frac{(2m+1)\pi}{2^{k+1}}$ for integer $m$, and $0\leq k\leq n$, with corresponding Laurent expansion \begin{equation} \begin{split} \mathcal{S}(x)&=6\cdot 16^k\left[\frac1{2^{4k}(x-x_0)^4}-\frac4{3\cdot 2^{2k}(x-x_0)^2}\right] -8\cdot 16^k\left[\frac1{2^{2k}(x-x_0)^2}\right]+O(1)\\ &=\frac{6}{(x-x_0)^4}+O(1) \end{split} \end{equation} about $x_0$. Using the same logic as above, we may construct the meromorphic, $\pi$-periodic function $$f(x)=16^{n+1}[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]-[6\cot^4(x)+8\cot^2(x)+2]$$ which also has a set of poles given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole $x_0$ also has the corresponding Laurent expansion $f(x)=\frac{6}{(x-x_0)^4}+O(1)$ about $x_0$.

Finally, notice that $\lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\tan^2(x)=-1$, so \begin{equation} \lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\mathcal{S}(x)-f(x)=\sum_{k=0}^n16^k[6-8+2]-16^{n+1}[6-8+2]+[6-8+2]=0 \end{equation}

Since $\mathcal{S}$ and $f$ are both meromorphic and $\pi$-periodic, and share the same poles with the same Laurent expansions about those poles up to $O(1)$, then $\mathcal{S}-f$ is analytic and $\pi$-periodic, and since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $S-f$ is bounded. Therefore by Louville's theorem, $\mathcal{S}-f$ is constant. Finally, since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $\mathcal{S}-f$ is precisely $0$, which completes the proof.

Addendum: Note that this procedure generalizes, and can be used to construct and prove higher order identities of this type. In fact, we have the following theorem which may be proven in the same way as above:

Let $c_{k,\ell}$ denote the $(-2k)$-th coefficient of the Laurent expansion of $\cot^{2\ell}(x)$ about $0$, and define the polynomial $P_m(u)=a_m u^{2m}+\cdots +a_1u^2+a_0$, where $a_0,a_1,\ldots, a_m$ is the solution to \begin{equation} \left[\begin{matrix} 0\\0\\0\\\vdots\\0\\1\\ \end{matrix}\right] = \left[\begin{matrix} 1&-1&1&\cdots&(-1)^{m-1}&(-1)^m\\ 0&c_{1,1}&c_{1,2}&\cdots&c_{1,m-1}&c_{1,m}\\ 0&0&c_{2,2}&\cdots&c_{2,m-1}&c_{2,m}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&c_{m-1,m-1}&c_{m-1,m}\\ 0&0&0&\cdots&0&c_{m,m}\\ \end{matrix}\right] \left[\begin{matrix} a_0\\a_1\\a_2\\\vdots\\a_{m-1}\\a_m\\ \end{matrix}\right] \end{equation} We then have the following identity: \begin{equation} \sum_{k=0}^n2^{2mk}P_m(\tan(2^k x))=2^{2m(n+1)}P_m(\cot(2^{n+1}x))-P_m(\cot(x)) \end{equation} For example, the first few polynomials $P_m(u)$ are given below: $$\begin{split} P_1(u)&=u^2+1\\ P_2(u)&=u^4+\frac{4}{3}u^2+\frac{1}{3}\\ P_3(u)&=u^6+2u^4+\frac{17}{15}u^2+\frac{2}{15}\\ P_4(u)&=u^8+\frac{8}{3}u^6+\frac{12}{5}u^4+\frac{248}{315}u^2+\frac{17}{315}\\ \end{split}$$ Note that the identity proven in the first part of this answer is a special case of the above more general identity, setting $m=2$ (just multiply both sides by $6$).

We shall prove the identity \begin{equation} \sum_{k=0}^n16^k[6\tan^4(2^kx)+8\tan^2(2^kx)+2]=[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]- 16^{n+1}[6\cot^4(x)+8\cot^2(x)+2] \end{equation} via complex analysis.

Let $\mathcal{S}=3S_1+2S_2+2S_3$ denote the LHS of the above equality, and notice that $\mathcal{S}$ is meromorphic, and $\pi$-periodic. Note that the poles of $\tan^4(x)$ and $\tan^2(x)$ lie at points of the form $x_0=(m+\frac12)\pi$, with corresponding Laurent expansions \begin{equation} \begin{split} \tan^4(x)&=\frac1{(x-x_0)^4}-\frac4{3(x-x_0)^2}+O(1)\\ \tan^2(x)&=\frac1{(x-x_0)^2}+O(1)\\ \end{split} \end{equation} about $x_0$. Therefore the set of poles of $\mathcal{S}$ is given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole may be written uniquely as $x_0=\frac{(2m+1)\pi}{2^{k+1}}$ for integer $m$, and $0\leq k\leq n$, with corresponding Laurent expansion \begin{equation} \begin{split} \mathcal{S}(x)&=6\cdot 16^k\left[\frac1{2^{4k}(x-x_0)^4}-\frac4{3\cdot 2^{2k}(x-x_0)^2}\right]-8\cdot 16^k\left[\frac1{2^{2k}(x-x_0)^2}\right]+O(1)\\ &=\frac{6}{(x-x_0)^4}+O(1) \end{split} \end{equation} about $x_0$. Using the same logic as above, we may construct the meromorphic, $\pi$-periodic function $$f(x)=[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]-16^{n+1}[6\cot^4(x)+8\cot^2(x)+2]$$ which also has a set of poles given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole $x_0$ also has the corresponding Laurent expansion $f(x)=\frac{6}{(x-x_0)^4}+O(1)$ about $x_0$.

Finally, notice that $\lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\tan^2(x)=-1$, so \begin{equation} \lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\mathcal{S}(x)-f(x)=\sum_{k=0}^n16^k[6-8+2]-[6-8+2]+16^{n+1}[6-8+2]=0 \end{equation}

Since $\mathcal{S}$ and $f$ are both meromorphic and $\pi$-periodic, and share the same poles with the same Laurent expansions about those poles up to $O(1)$, then $\mathcal{S}-f$ is analytic and $\pi$-periodic, and since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $S-f$ is bounded. Therefore by Louville's theorem, $\mathcal{S}-f$ is constant. Finally, since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $\mathcal{S}-f$ is precisely $0$, which completes the proof.

We shall prove the identity \begin{equation} \sum_{k=0}^n16^k[6\tan^4(2^kx)+8\tan^2(2^kx)+2]=16^{n+1}[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]- [6\cot^4(x)+8\cot^2(x)+2] \end{equation} via complex analysis.

Let $\mathcal{S}=3S_1+2S_2+2S_3$ denote the LHS of the above equality, and notice that $\mathcal{S}$ is meromorphic, and $\pi$-periodic. Note that the poles of $\tan^4(x)$ and $\tan^2(x)$ lie at points of the form $x_0=(m+\frac12)\pi$, with corresponding Laurent expansions \begin{equation} \begin{split} \tan^4(x)&=\frac1{(x-x_0)^4}-\frac4{3(x-x_0)^2}+O(1)\\ \tan^2(x)&=\frac1{(x-x_0)^2}+O(1)\\ \end{split} \end{equation} about $x_0$. Therefore the set of poles of $\mathcal{S}$ is given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole may be written uniquely as $x_0=\frac{(2m+1)\pi}{2^{k+1}}$ for integer $m$, and $0\leq k\leq n$, with corresponding Laurent expansion \begin{equation} \begin{split} \mathcal{S}(x)&=6\cdot 16^k\left[\frac1{2^{4k}(x-x_0)^4}-\frac4{3\cdot 2^{2k}(x-x_0)^2}\right] -8\cdot 16^k\left[\frac1{2^{2k}(x-x_0)^2}\right]+O(1)\\ &=\frac{6}{(x-x_0)^4}+O(1) \end{split} \end{equation} about $x_0$. Using the same logic as above, we may construct the meromorphic, $\pi$-periodic function $$f(x)=16^{n+1}[6\cot^4(2^{n+1}x)+8\cot^2(2^{n+1}x)+2]-[6\cot^4(x)+8\cot^2(x)+2]$$ which also has a set of poles given by $(\frac{\pi}{2^{n+1}})\mathbb{Z}\backslash \pi\mathbb{Z}$, and each pole $x_0$ also has the corresponding Laurent expansion $f(x)=\frac{6}{(x-x_0)^4}+O(1)$ about $x_0$.

Finally, notice that $\lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\tan^2(x)=-1$, so \begin{equation} \lim_{\mathfrak{I}(x)\rightarrow\pm\infty}\mathcal{S}(x)-f(x)=\sum_{k=0}^n16^k[6-8+2]-16^{n+1}[6-8+2]+[6-8+2]=0 \end{equation}

Since $\mathcal{S}$ and $f$ are both meromorphic and $\pi$-periodic, and share the same poles with the same Laurent expansions about those poles up to $O(1)$, then $\mathcal{S}-f$ is analytic and $\pi$-periodic, and since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $S-f$ is bounded. Therefore by Louville's theorem, $\mathcal{S}-f$ is constant. Finally, since $\mathcal{S}-f$ vanishes as $\mathfrak{I}(x)\rightarrow\pm\infty$, then $\mathcal{S}-f$ is precisely $0$, which completes the proof.

Addendum: Note that this procedure generalises, and can be used to construct and prove higher order identities of this type. In fact, we have the following theorem which may be proven in the same way as above:

Let $c_{k,\ell}$ denote the $(-2k)$-th coefficient of the Laurent expansion of $\cot^{2\ell}(x)$ about $0$. If $a_1,\ldots, a_m$ is the solution to \begin{equation} \left[\begin{matrix} 0\\0\\0\\\vdots\\0\\1\\ \end{matrix}\right] = \left[\begin{matrix} 1&-1&1&\cdots&(-1)^{m-1}&(-1)^m\\ 0&c_{1,1}&c_{1,2}&\cdots&c_{1,m-1}&c_{1,m}\\ 0&0&c_{2,2}&\cdots&c_{2,m-1}&c_{2,m}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\cdots&c_{m-1,m-1}&c_{m-1,m}\\ 0&0&0&\cdots&0&c_{m,m}\\ \end{matrix}\right] \left[\begin{matrix} a_0\\a_1\\a_2\\\vdots\\a_{m-1}\\a_m\\ \end{matrix}\right] \end{equation} then \begin{equation} \sum_{k=0}^n2^{2mk}\sum_{\ell=0}^ma_\ell\tan^{2\ell}(2^k x)=\sum_{\ell=0}^m a_\ell[2^{2m(n+1)}\cot^{2\ell}(2^{n-1}x)-\cot^{2\ell}(x)] \end{equation}