We can derive the normalizing constant $C(\lambda):=\int e^{-\lambda ||M||_*}dM$ using the smooth coarea formula.
Suppose we are given smooth functions $F: \mathbb{R}^n\to\mathbb{R}$ and $f: \mathbb{R}\to\mathbb{R}$. The formula states that
$$\int_{\mathbb{R}^n} f(F(x))dx=\int_{t\in \mathbb{R}}\left(\int_{x\in F^{-1}t} f(x) |\nabla F(x)|^{-1}d H_t(x) \right)dt$$
where $dH_t$ denotes the surface measure on the preimage $F^{-1}t$. In our case, we'll have $f(x)=e^{-\lambda x}$ and $F(M)=||M||_*$. Also, denote the dimensions of $M$ by $n,m$
So first of all we'll need to compute $|\nabla F(M)|$. By the resuls from these notes, we can write down the derivatives of the singular values of a matrix. Suppose we have an SVD decomposition $M=USV^t$, with $U^tU=I_{min(n,m)}=V^tV$ and $S$ the diagonal matrix of singular values (we can assume $M$ is full rank, because the corresponding set of matrices has full measure).Then the differential with respect to the ith singular value is given by
$dS_i=(U^tdMV)_{ii}$. Summing over these differentials, and recalling that the nuclear norm is the sum of singular values, we get
\begin{eqnarray*}
d||M||_{*}&=&Tr(U^tdMV)\\
&=&\sum_i \sum_{a,b}U[a,i]V[b,i]dM_{ab}\\
& = & \sum_{a,b} dM_{ab} \sum_i U[a,i]V[b,i]\\
& = & \sum_{ab}dM_{ab} (UV^t)_{ab}
\end{eqnarray*}
In other words, if we consider the gradient $\nabla ||M||_*$ with respect to the entries of $M$, then the $ab$ coefficient is given by $(UV^t)_{ab}$. Now computing the norm of the gradient is straight forward:
\begin{eqnarray*}
|\nabla ||M||_*|^2& = & \sum_{ab} (UV^t)_{ab}^2\\
& = & Tr((UV^t)(UV^t)^t)\\
& = & Tr(UV^tVU^t)\\
& = & Tr(U^tUV^tV)\\
& = & Tr(I_{min(n,m)}I_{min(n,m)})\\
& = & min(n,m)
\end{eqnarray*}
In the second equality we used the fact that $\sum_{ab} X_{ab}^2=Tr(XX^t)$ for any matrix $X$.
Now, we can compute $\int e^{-\lambda F(M)}dM$ using the coarea formula (and recall that $F(M):=||M||_*$). We get
$$\int_0^{\infty} \left(\int_{M\in B^*(t)} |\nabla F(M)|^{-1}e^{-\lambda F(M)} dM_{|B^*(t)}\right)dt$$
where $B^*(t)$ is the set of all matrices of nuclear norm $t$, and $dM_{|B^*(t)}$ denotes the restriction of the Euclidean measure to this set. Conveniently, the integrand is constant over each preimage, $$min(n,m)^{-1/2}\int_0^{\infty} e^{-\lambda t} \int_{B^*(t)} dM_{|B^*(t)}dt=min(n,m)^{-1/2}\int_0^{\infty} e^{-\lambda t} Vol(B^*(t))$$
The ball $B^*(t)$ is a dilation of the unit ball $B^*(1)$,so the volume is $Vol(B^*(1)) t^{nm-1}$.
Finally, we get
$$min(n,m)^{-1/2}Vol(B^*(1))\int_0^{\infty}e^{-\lambda t}t^{nm-1}dt=min(n,m)^{-1/2}Vol(B^*(1))(nm-1)!\lambda^{-nm}$$
for the normalizing constant $C(\lambda)$.
Given this, we can also derive the distribution of the random variable $||M||_*$. To start with, we'll take the identity $\int e^{-\lambda ||M||_*}=C(\lambda)$ and differentiate under the integral to see:
$(-1 )^k\int ||M||^k_* e^{-\lambda ||M||_*}=d^k_{\lambda}C(\lambda)=(-1)^k{\frac {(nm)^{(k)}}{\lambda^k}}C(\lambda)$
where $(nm)^{(k)}={\frac {(nm+k-1)!}{(nm-1)!}}$ denotes the rising factorial. Thus
$$E ||M||^k_*=C(\lambda)^{-1}\int ||M||^k_* e^{-\lambda ||M||_*}dM= \lambda^{-k}{\frac {(nm+k-1)!}{(nm-1)!}}=\lambda^{-k}{\frac {\Gamma(nm+k)}{\Gamma(nm)}}$$
Remarkably, these are the same moments as a $\Gamma$ distribution with rate parameter $\lambda$ and shape parameter $mn$. By uniqueness of moment generating functions, we conclude that $||M||_*\sim \Gamma(mn,\lambda)$.
added later
You can also use a similar argument to derive the distribution over singular vectors/values. To phrase the question, suppose that I have sampled orthogonal matrices $U,V$ independently and uniformly at random, as well as a positive vector $\sigma$ from some fixed distribution $p(x)$ on $\mathbb{R}_+^n$. The question is: What should $p$ be, in order to ensure that the resulting matrix $M:=UDiag(\sigma)V^t$ is distributed according to the original distribution $e^{-\lambda ||M||_*}$?
Assume WLOG that $m>n$. By the change of variables formula, the resulting density of $M$ is given by $|{\frac {\partial (U,V,\sigma)}{\partial M}}|p(\sigma)$,were $U,V,\sigma$ are the SVD of $M$.
Lemma
$$\left|{\frac {\partial (U,V,\sigma)}{\partial M}}\right|=\prod_{i<=n} \sigma_i^{-(m-n)}\prod_{i<j\leq n}{\frac 1 {|\sigma_i^2-\sigma_j^2|}}$$
proof
See Aspects of Multivariate Statistical Theory by R.J. Muirhead. I also worked it out in a previous edit to this question, so it is visible in the edit history. It is not particularly difficult, provided that one uses the formulas for the derivative of the SVD in the linked notes. There are two tricks that make it considerably easier: (1) by symmetry it suffices to assume that $U=I$ and $V=\left (\begin{array}{c}I\\ 0\end{array}\right)$, (2) rather than the raw jacobian, it turns out to be easier to work with the Gram matrix $G_{ab}=\sum_{ij} {\frac {\partial x_a}{\partial M_{ij}}}{\frac {\partial x_b}{\partial M_{ij}}}$ where $x_a$ ranges over all independent entries of $U,V,\sigma$. This matrix ends up having a particularly simple structure.
So in summary:
The original distribution $e^{-\lambda ||M||_*}$ is obtained as the distribution of a matrix of the form $UDiag(\sigma)V^T$, where $U$ and $V$ are sampled uniformly from the space of all orthogonal matrices, and $\sigma$ is sampled from $\mathbb{R}_+^n$ according to the density $\propto \prod_{i\leq n} \sigma_i^{m-n}\prod_{i<j\leq n}|\sigma_i^2-\sigma_j^2| e^{-\lambda \sum_i\sigma_i}$.
singular value distribution
To explore the singular values, I generated samples from the singular value distribution using MCMC. in this case, $n=m=20$ and $\lambda=1$. Below is the histogram of the resulting sums $\sum_i \sigma_i$. These sums appear to follow a $\Gamma(20*20,1)$ distribution (shown in orange), as predicted by the above discussion.
A bit more interesting is the distribution over the individual singular values. If $\sigma_k$ denotes the $k$th smallest singular value, then it appears that each $\sigma_k$ is in fact well approximated by a gamma distribution. Below are plotted histograms for a range of $\sigma_k$, with the approximating gamma pdfs overlaid.
This is quite interesting, and I am not sure why it would be the case.
analysis of largest singular value
It turns out that there is a relatively simple analytic expression for the largest singular value, which can be found by applying the techniques from this paper. The key property is that the singular value density looks like $\propto |det \{\phi_j(\sigma_i)\}_{ij}|$ where $\phi_i$ are scalar-valued functions. Due to the Vandermonde determinant formula, it is easy to see that we can take $\phi_j(\sigma)=\sigma^{m-n}e^{-\lambda\sigma}\sigma^{2(j-1)}$.
To evaluate the cdf of the largest singular value, we need to evaluate the integral $\int_{0<\sigma_1<\dots<\sigma_n<t}|det \{\phi_j(\sigma_i)\}_{ij}|d\sigma$. In general, it turns out that an integral of this form can be expressed as $\sqrt{|det(A)|}$ where $A_{ij}=\int_{[0,t]^2}sgn(x-y)\phi_i(x)\phi_j(t)dxdy$ (for even $n$; if n is odd the definition of $A$ is slightly more complicated). Even better, it turns out that, for our $\phi$ functions, the entries of $A$ satisfy a recurrence relation that obviates the need to actually evaluate the integral. See the linked paper for more details.
