  In this work, we investigated Bayesian methodologies for constraining in the
Solar System a Yukawa suppression of the Newtonian potential -- which we
interpret as the effect of a non-null graviton mass -- by considering its
impact on planetary orbits. Complementary to the previous results obtained with
INPOP planetary ephemerides, we consider here a Markov Chain Monte Carlo
approach associated with a Gaussian Process Regression for improving the
resolution of the constraints driven by planetary ephemerides on the graviton
mass in the Solar System. At the end of the procedure, a posterior for the mass
of the graviton is presented, providing an upper bound at $1.01 \times 10^{-24}
\; eV c^{-2}$ (resp. $\lambda_g \geq 122.48 \times 10^{13} \; km$) with a
$99.7\%$ confidence level. The threshold value represents an improvement of 1
order of magnitude relative to the previous estimations. This updated
determination of the upper bound is mainly due to the Bayesian methodology,
although the use of new planetary ephemerides (INPOP21a used here versus
INPOP19a used previously) already induces a gain of a factor 3 with respect to
the previous limit. The INPOP21a ephemerides is characterized by the addition
of new Juno and Mars orbiter data, but also by a better Solar System modeling,
with notably a more realistic model of the Kuiper belt. Finally, by testing the
sensitivity of our results to the choice of the $\textit{a priori}$
distribution of the graviton mass, it turns out that the selection of a prior
more favorable to zero-mass graviton (that is, here, General Relativity) seems
to be more supported by the observations than non-zero mass graviton, leading
to a possible conclusion that planetary ephemerides are more likely to favor
General Relativity.
