I can't justify which method is correct and why. If I have an $(M x N)$ matrix and I want to calculate its covariance, would mean center the matrix and then do
$$ X^TX $$
or...
$$ XX^T $$
does the choice of which one to do depend on whether I want to get a covariance of the rows or the columns?
It all depends on what the rows and the columns represent. If your columns correspond to different dimensions/predictors and each row to different instances, then the correlation for 2 dimension/predictor pairs ($n_{1}$,$n_{2}$) is just
\begin{align} \operatorname{cov}(X[,n_1], X[,n_2]) &= \operatorname{E}\left[X[,n_1] X[,n_2]\right] - \operatorname{E}\left[X[, n_1]\right] \operatorname{E}\left[X[,n_2]\right], \\ &= \operatorname{E}\left[X[,n_1] X[,n_2]\right] \\ \end{align}
Where in the last equality the means of columns $n_1$ and $n_2$ are taken to be $0$. Thus, in this case the correlation matrix would be $X^{T}X$. In the case where your columns represent instances and your rows different dimensions/predictors, then $XX^{T}$.