This is a complement to my other answer to this question.
The other answer gives an explicit basis when $n=3$.
This answer shows the generation property (but does not give an explicit basis) when $n\geq 3$.
I thought it best to post it as a separate answer.
For any $n\geq 3$, ${\cal M}_{n,n}(\mathbb R)$ is indeed generated by
$SO(n)$, it is in fact already generated by the signed permutation matrices
whose determinant is $+1$ (let us denote by $Z(n)$ the set of all such
matrices).
(Recall that a signed permutation matrix is a matrix which has exactly
one nonzero coefficient in each row and column, and this coefficient
is always $\pm 1$).
To check this, denote by $W$ the subspace generated by those matrices.
Denote by $E_{ij}$ the matrix all of those coefficents are zero except
the one in place $(i,j)$, equal to 1 (so $(E_{ij})_{1\leq i,j \leq n}$
is the canonical basis of ${\cal M}_{n,n}(\mathbb R)$).
Any two $E_{i,i}$ are conjugate, and any two $E_{i,j} (i\neq j)$ are
conjugate. Since $SO(n)$ (and hence $W$ also) are invariant by
conjugation, it will suffice to show that $E_{1,1}$ and $E_{1,2}$ are
in $W$.
In fact, it will suffice to show that $E_{1,2}\in W$, because once
we know that $E_{1,2}\in W$, we have an element
$T=(t_{ij})\in Z(n)$ with $t_{11}=1,t_{j(j+1)}=\pm 1$
(for $2\leq j \leq n-1$),$t_{n2}=\pm 1$. Then $E_{1,1}-T$ is a sum
of precisely $n-1$ matrices of the form $E_{xy}$ with $x\neq y$ ; and we deduce
that $E_{1,1} \in T$ also.
So let us show that $E_{1,2}\in W$.
For any sequence of signs $\varepsilon=
(\varepsilon_1,\varepsilon_2, \ldots,
\varepsilon_n) \in \lbrace\pm 1\rbrace ^n$, consider the following two
signed permutation matrices :
$$
A(\varepsilon)=\left(
\begin{array}{ccccccc}
0 & \varepsilon_1 & 0 & 0 & \ldots & 0 & 0 \\
0 & 0 & \varepsilon_2 & 0 & \ldots & 0 & 0 \\
0 & 0 & 0 & \varepsilon_3 & \ldots & 0 & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & 0 & \ldots & \varepsilon_{n-2} & 0\\
0 & 0 & 0 & 0 & \ldots & 0 & \varepsilon_{n-1}\\
\varepsilon_n & 0 & 0 & 0 & \ldots & 0 & 0 \\
\end{array}\right) ,
B(\varepsilon)=\left(
\begin{array}{ccccccc}
0 & \varepsilon_1 & 0 & 0 & \ldots & 0 & 0 \\
0 & 0 & \varepsilon_2 & 0 & \ldots & 0 & 0 \\
0 & 0 & 0 & \varepsilon_3 & \ldots & 0 & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & 0 & \ldots & \varepsilon_{n-2} & 0\\
\varepsilon_n & 0 & 0 & 0 & \ldots & 0 & 0 \\
0 & 0 & 0 & 0 & \ldots & 0 & \varepsilon_{n-1}\\
\end{array}\right)
$$
so that $B(\varepsilon)$ has been obtained from $A(\varepsilon)$ by swapping
the last two rows. If we put $p=\prod_{k}\varepsilon_k$ , the determinant
of $A(\varepsilon)$ is $(-1)^{n-1}p$, and the determinant of $B$ is
$(-1)^np$.
When $n$ is odd, we simply have
$E_{12}=\frac{B(1,-1,-1,\ldots,-1)+B(1,1,1,\ldots,1)}{2}$.
When $n$ is even, things are slightly more complicated : if we put
$$C=\bigg\lbrace (\varepsilon_2,\varepsilon_3, \ldots,
\varepsilon_{n}) \in \lbrace\pm 1\rbrace ^{n-1}
\bigg| \ \prod_{k=2}^{n}\varepsilon_k=(-1)\bigg\rbrace,
$$
$$
C'=\bigg\lbrace( 1,\varepsilon_2,\varepsilon_3, \ldots,
\varepsilon_{n}) \bigg| (\varepsilon_2,\varepsilon_3, \ldots,
\varepsilon_{n}) \in C \bigg\rbrace
$$
then
$$
E_{12}=\frac{\displaystyle\sum_{\varepsilon \in C'}A(\varepsilon)}{2^{n-1}}
$$
In both cases, we have written $E_{12}$ as a linear combination
of matrices in $Z(n)$, so we are done.