I don't have enough reputation to comment, but for what it's worth, if we replace $\tan(\theta_i)$ with $\lambda_i$ and set $P(x)=\prod (x-\lambda_i)$, then the left hand side will be $P(i)P(-i)$ by difference of squares and Vieta's formula, and the right hand side $\prod (1+\lambda_i^2)$. The equality then follows from $(i-\lambda_i)(-i-\lambda_i)=1+\lambda_i^2$. In this case, all that was used was $\tan^2(\theta)+1=\sec^2(\theta)$.

I don't have enough reputation to comment, but for what it's worth, if we replace $\tan(\theta_j)$ with $\lambda_j$ and set $P(x)=\prod (x-\lambda_j)$, then the left hand side will be $P(i)P(-i)$ by difference of squares and Vieta's formula, and the right hand side $\prod (1+\lambda_j^2)$. The equality then follows from $(i-\lambda_j)(-i-\lambda_j)=1+\lambda_j^2$. In this case, all that was used was $\tan^2(\theta)+1=\sec^2(\theta)$.