Ursell functions $U_n$ are higher-order generalizations of the covariance function, which capture the interactions between $n$ random variables. In the classical Ising model, as shown by Shlosman, when considering the spins at some locations, the sign of $U_{2n}$ alternates with $n$ and is independent of the locations of the spins considered. In this paper, we study the Ursell function in Ising lattice gauge theory. When the spins at the edges are used as random variables, we show that $U_n$ can be positive, negative, or zero depending on the configuration and the parameter $β$. When considering Wilson loops observables as random variables, using the tool of cluster expansion adapted to this setting, we prove that at sufficiently low temperature, for any number $n$ of disjoint Wilson loops, there exists a configuration of loops such that the Ursell function $U_n$ is positive. These results contrast sharply with the behavior observed for the Ising model.