I'm guessing that you meant to write "real orthogonal matrices" without "symmetric", since the set of symmetric orthogonal matrices has Haar measure 0. Otherwise, please clarify what you mean by "a distribution consistent with Haar measure over $O(N)$".
Individual entries of Haar-distributed $O(N)$ matrices are approximately Gaussian; in some form this goes back to Maxwell. This has recently been improved in two slightly different directions. First, fairly large submatrices are approximately distributed like matrices with i.i.d. Gaussian entries; see this paper by Jiang. Second, linear combinations of the entries are approximately Gaussian, in a very strong sense; see this paper by E. Meckes. They are even approximately jointly Gaussian; see this paper by Chatterjee and E. Meckes.
Update: I discussed this question with Elizabeth, who pointed out to me that my suggestion isn't necessarily great because Gaussian behavior can creep in for all sorts of reasons, in particular due to an accumulation of small, nearly independent errors (i.e., the good old central limit theorem). But you the suggestion may be salvageable by getting more quantitative. For example, Elizabeth's paper shows that linear functionals approximate Gaussians with an error (in total variation) of $c/N$, which is better than one would expect from CLT effects. Even better, you can consider traces of powers, which were shown by Johansson to converge exponentially quicky to a Gaussian distribution.
