We introduce a distributionally robust maximum likelihood estimation model\nwith a Wasserstein ambiguity set to infer the inverse covariance matrix of a\n$p$-dimensional Gaussian random vector from $n$ independent samples. The\nproposed model minimizes the worst case (maximum) of Stein's loss across all\nnormal reference distributions within a prescribed Wasserstein distance from\nthe normal distribution characterized by the sample mean and the sample\ncovariance matrix. We prove that this estimation problem is equivalent to a\nsemidefinite program that is tractable in theory but beyond the reach of\ngeneral purpose solvers for practically relevant problem dimensions $p$. In the\nabsence of any prior structural information, the estimation problem has an\nanalytical solution that is naturally interpreted as a nonlinear shrinkage\nestimator. Besides being invertible and well-conditioned even for $p>n$, the\nnew shrinkage estimator is rotation-equivariant and preserves the order of the\neigenvalues of the sample covariance matrix. These desirable properties are not\nimposed ad hoc but emerge naturally from the underlying distributionally robust\noptimization model. Finally, we develop a sequential quadratic approximation\nalgorithm for efficiently solving the general estimation problem subject to\nconditional independence constraints typically encountered in Gaussian\ngraphical models.\n