Closed, sum-free form for the $n$-th derivative of $\operatorname{arcsinh}(\frac1x)$ in $x=1$
Here is a proof that $M_{nn}=O(n^{-1})$, as the OP conjectured. The function $$f(z):=\operatorname{arcsinh}(1/z),\qquad |z|<1,$$ is holomorphic in the unit disk, hence by Cauchy's formula $$ \frac{f^{(k)}(1)}{k!}=\frac{1}{2\pi i}\int_{|z-1|=r}\frac{f(z)}{(z-1)^{k+1}}\,dz,\qquad 0<r<1.$$ Here and later circles are positively oriented. It follows, by applying the binomial formula, that $$\frac{1}{2}M_{nn}=\frac{1}{2\pi i}\int_{|z-1|=r}\frac{f(z)}{z-1}\left(\frac{z+1}{z-1}\right)^n\,dz,\qquad 0<r<1.$$ We make the change of variables $z:=(1+w)/(1-w)$, keep track of the image of the original circle (which is another circle), and deform it to a circle centered at the origin. This way we obtain $$\frac{1}{2}M_{nn}=\frac{1}{2\pi i}\int_{|w|=\rho}\frac{g(w)}{1-w}\cdot\frac{dw}{w^{n+1}},\qquad 0<\rho<1,$$ where $g(w)$ is the following holomorphic function in the unit disk: $$g(w):=\operatorname{arcsinh}\left(\frac{1-w}{1+w}\right),\qquad |w|<1.$$ This way we see, again by Cauchy's formula, that $M_{nn}$ is twice the sum of the first $n$ Taylor coefficients of $g(w)$ at the origin. (Added: This can also be derived without complex analysis, as in Max Alekseyev's response.) Using the identities $$g'(w)=\frac{-\sqrt{2}}{(1+w)(1+w^2)^{1/2}},\qquad|w|<1,$$ $$\frac{1}{(1+w^2)^{1/2}}=\sum_{k=0}^\infty\frac{(-1)^k}{4^k}\binom{2k}{k}w^{2k},\qquad|w|<1,$$ it is straightforward to derive that $$\frac{1}{2}M_{nn}=\operatorname{arcsinh}(1)-\sqrt{2}\sum_{2k\leq n}\frac{(-1)^k}{4^k}\binom{2k}{k}\sum_{2k\leq m\leq n-1}\frac{(-1)^m}{m+1}.$$ The inner sum is positive and it is of size $O(k^{-1})$, hence in the outer sum the $k$-th term has sign $(-1)^k$ and it is of size $O(k^{-3/2})$. It follows that $M_{nn}$ converges. By applying a more careful asymptotic analysis in $k$ and keeping track of the error term $O(n^{-1})$ coming from the inner sum, the convergence is seen to occur with speed $O(n^{-1})$. Finally, by Abel's theorem, $$\lim_{n\to\infty} M_{nn}=2g(1)=0,\qquad\text{hence in fact}\qquad M_{nn}=O(n^{-1}).$$
The numbers $c_k = \frac{1}{k!} \left.\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x\right|_{x=1}$ are the coefficients in the expansion: $$\operatorname{arcsinh}\frac1x = c_0 + c_1(x-1) + c_2(x-1)^2 + \dots.$$ It follows that $2^kc_k$ is the coefficient of $t^k$ in $\operatorname{arcsinh}\frac1{1+2t}$. Now, $$M_{nn} = 2\sum_{k=0}^n \binom{n}{k} 2^kc_k = 2\cdot [t^n]\ (1+t)^n \operatorname{arcsinh}\frac1{1+2t}.$$ Using Lagrange inversion, we get the generating function for $M_{nn}$: $$\sum_{n\geq 0} M_{nn} t^n = \frac{2}{1-t} \operatorname{arcsinh}\frac{1-t}{1+t}.$$ From here the asymptotic for $M_{nn}$ can be obtained using the standard tools (e.g., see the answer from GH from MO).
To get an explicit formula for $\left.\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x\right|_{x=1}$ (and thus for $c_k$), we notice that $$\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x = - \left(\frac{d}{dx}\right)^{k-1} (x^2+x^4)^{-\frac12}.$$ To expand the last expression, one can use Faà di Bruno's formula: $$\left(\frac{d}{dx}\right)^{k-1} (x^2+x^4)^{-\frac12}$$ $$ = \sum_{i=1}^{k-1} \binom{-\frac{1}{2}}{i} i! (x^2+x^4)^{-\frac12-i} B_{k-1,i}(2x+4x^3,2+12x^2,24x,24,0,0,\dots,0),$$ where $B_{k-1,j}$ are Bell polynomials.
Evaluating at $x=1$, for $k>0$, we get $$\left.\left(\frac{d}{dx}\right)^k\operatorname{arcsinh}\frac1x\right|_{x=1} = -\sum_{i=1}^{k-1} \binom{-\frac{1}{2}}{i} i! 2^{-\frac12-i} B_{k-1,i}(6,14,24,24,0,0,\dots,0)$$ $$ = -\frac{(k-1)!}{\sqrt{2}}\sum_{j_1+2j_2+3j_3+4j_4=k-1} \frac{(-1)^{j_1+j_2+j_3+j_4}(2(j_1+j_2+j_3+j_4))!}{(j_1+j_2+j_3+j_4)!j_1!j_2!j_3!j_4!} 2^{-3(j_1+j_2+j_3+j_4)} 6^{j_1}7^{j_2}4^{j_3}$$ $$=-\frac{(k-1)!}{\sqrt{2}}\sum_{j_1+2j_2+3j_3+4j_4=k-1} \frac{(-1)^{j_1+j_2+j_3+j_4} (2(j_1+j_2+j_3+j_4))!}{(j_1+j_2+j_3+j_4)!j_1!j_2!j_3!j_4!} 2^{-2j_1-3j_2-j_3-3j_4} 3^{j_1}7^{j_2}.$$