We consider the additive version of the matrix denoising problem, where a random symmetric matrix S of size n has to be inferred from the observation of \(Y=S+Z\) , with Z an independent random matrix modeling a noise. For prior distributions of S and Z that are invariant under conjugation by orthogonal matrices we determine, using results from first and second order free probability theory, the Bayes-optimal (in terms of the mean square error) polynomial estimators of degree at most D, asymptotically in n, and show that as D increases they converge towards the estimator introduced by Bun et al. (IEEE Trans Inf Theory 62:7475, 2016). We conjecture that this optimality holds beyond strictly orthogonally invariant priors, and provide partial evidences of this universality phenomenon when S is an arbitrary Wishart matrix and Z is drawn from the Gaussian Orthogonal Ensemble, a case motivated by the related extensive rank matrix factorization problem.