Discretize the time domain and regard the time-sequences as $N\times 1$ matrices; then you have $\mathbf s_1^T~\mathbf s_2=0.$ This product is preserved by any invertible linear transform $\mathbf s\mapsto \mathbf C~\mathbf s$ as long as $\mathbf C^T~\mathbf C = \mathbf I,$ the orthogonal group. Complexify this by converting $\bullet^T\mapsto\bullet^\dagger$and you get a nice result; the Fourier transform looks like,
$$s[f] =\mathcal F_{t\to f}~s(t) = \int_{-\infty}^\infty dt~e^{-2\pi ift}~s(t),\\s(t) =\mathcal F^{-1}_{f\to t}~s[f] = \int_{-\infty}^\infty df~e^{2\pi ift}~s[f].$$
Note that under the transpose $f\leftrightarrow t$ you don't actually see this $\mathcal F$ matrix become its inverse (it remains identical; it's a symmetric matrix), but upon conjugation $i \leftrightarrow -i$ it definitely does. So your Fourier matrix is a unitary matrix $\mathbf F^\dagger~\mathbf F = 1.$ (Actually looking over the discrete Fourier transform formula, you may need to have a constant $\alpha$ between them: but it looks like it's just a real number and so you can divide through and normalize appropriately.)
You can therefore see that $(\mathbf F ~ \mathbf s_1)^\dagger~(\mathbf F ~ \mathbf s_2) = \mathbf s_1^\dagger~\mathbf s_2 = 0,$ if and only if the two things are time-orthogonal. The cost is that you need to instead think of frequency-orthogonality as being $$\langle s_1,~s_2\rangle_f = \int df~s_1^*[f]~s_2[f].$$If you're not willing to do that then the best I can give you is this: $\mathbf C^T \mathbf C$ is a symmetric matrix and therefore in any basis the notion of orthogonality in some other basis becomes some symmetric kernel $g(t, t')$ with that inner product coming into the time basis as$$\langle s_1, s_2\rangle_g = \iint dt~ dt'~g(t, t')~s_1(t)~s_2(t').$$Of course multiplication in the frequency domain is going to look like convolution in the time domain, so probably for the Fourier transform, $g$ is some sort of Dirac $\delta$-function.