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.