Assume that X and Y are topological spaces with $\sigma$-finite measure, $L_p(X)$ is a Banach space of complex-valued functions so that $\int_X |f(x)|^p dx < \infty$, $1 \leqslant p < \infty$. Assume that the set $\{e_1(x), e_2(x), ... e_i(x), ...\}_{i = 1}^\infty$ is an unconditional basis in $L_p(X)$, $\{\eta_1(x), \eta_2(x), ... \eta_j(x), ...\}_{j = 1}^\infty$ is an unconditional basis in $L_p(Y)$. Then the set $\{e_1 \eta_1, e_2 \eta_2, ..., e_1 \eta_j, ..., e_2 \eta_1, e_2, \eta_2, ..., e_2 \eta_j, ..., e_i \eta_1, e_i \eta_2, ..., e_i \eta_j, ...\}_{i, j = 1}^\infty$ is the unconditional basis in the $L_p(X \times Y)$ space.
It seems a well-known fact, but I failed to find where it is proven. It would be great if someone can help me to find a reference to this fact.