The group delay $\tau(\omega)$ of a Bessel filter equals $1/\omega_0$ at $\omega=0$, where $\omega_0$ is the normalization frequency. It is independent of the filter order, unless we choose different normalization frequencies for different filter orders (see below for the normalization used by Matlab and scipy).
Note that the value $\tau(0)=1/\omega_0$ is not an approximation but it is exact. Since Bessel filters are characterized by a maximally flat group delay at DC, the value $1/\omega_0$ is approximated in a maximally flat way throughout the passband. So for fixed $\omega_0$, the only difference between the group delays of Bessel filters with different orders is where exactly and how the group delay starts deviating from the value $1/\omega_0$. But that happens somewhere in the stopband, where the group delay is not relevant anyway.
The group delay of a Bessel filter can be expressed in terms of reverse Bessel polynomials, but the extra information one obtains is mainly about the group delay response in the stopband, because the passband group delay is flat and approximates $1/\omega_0$ very closely, as mentioned above.
The group delay of an $n^{th}$ order Bessel filter is given by
$$\tau_n(\omega)=\frac{1}{\omega_0}\textrm{Im}\left\{\frac{\theta_n'\left(j\frac{\omega}{\omega_0}\right)}{\theta_n\left(j\frac{\omega}{\omega_0}\right)}\right\}\tag{1}$$
$\theta_n(s)$ is an $n^{th}$ order reverse Bessel polynomial as specified here:
$$\theta_n(s)=\sum_{k=0}^{n}a_ks^k,\qquad a_k=\frac{(2n-k)!}{2^{n-k}k!(n-k)!}\tag{2}$$
$\theta_n'(j\Omega)$ is the derivative with respect to $\Omega$.
For $\omega=0$ we obtain from $(1)$
\begin{align*} \tau_n(0) &= \frac{1}{\omega_0}\textrm{Im}\left\{\frac{\theta_n'\left(0\right)}{\theta_n\left(0\right)}\right\}\\ &= \frac{1}{\omega_0}\textrm{Im}\left\{\frac{ja_1}{a_0}\right\}\\&= \frac{1}{\omega_0}\frac{a_1}{a_0}\\&= \frac{1}{\omega_0}\tag{3} \end{align*}
because from the definition of the coefficients $a_k$ in $(2)$ it follows that $a_0=a_1$ always holds, regardless of the order $n$.
Eq. $(3)$ shows that the group delay in the passband only depends on the chosen normalization frequency $\omega_0$. Note that the group delay is maximally flat at $\omega=0$, so it remains almost constant for at least a large portion of the passband, and - depending on the filter order - this almost constant behavior can even extend into the transition band or a part of the stopband.
Most implementations of design routines for Bessel filters use different normalization frequencies $\omega_0$ for different filter orders, so it appears that the group delay depends on the filter order, as observed by the OP. The normalization used in Matlab's besselap makes sure that the first and the last coefficients of the filter's denominator polynomial equal $1$. This can only be achieved if different normalization frequencies $\omega_0$ are chosen for different filter orders. The transfer function of an $n^{th}$ order Bessel filter is given by
\begin{align*} H(s) &= \frac{\theta_n(0)}{\theta_n\left(\frac{s}{\omega_0}\right)} \\ &= \frac{a_0}{\left(\frac{s}{\omega_0}\right)^n+\ldots+a_0}\\ &= \frac{1}{\frac{1}{a_0}\left(\frac{s}{\omega_0}\right)^n+\ldots+1}\tag{4} \end{align*}
From $(4)$, the leading coefficient of the denominator is given by
$$\frac{1}{a_0\omega_0^n}\tag{5}$$
It equals $1$ if
$$\omega_0=a_0^{-1/n}\tag{6}$$
Consequently, with this normalization, the normalization frequency $\omega_0$ depends on the filter order $n$, and the group delay at $\omega=0$ is
$$\tau_n(0)=\frac{1}{\omega_0}=a_0^{1/n}=\left(\frac{(2n)!}{2^nn!}\right)^{1/n}\tag{7}$$
which is consistent with what is claimed here.
According to the scipy-documentation, this is also the default normalization used in scipy.signal.bessel
(norm=phase
).
Summarizing, the group delay at DC of a Bessel filter equals the inverse normalization frequency. It is a common choice to use different normalization frequencies for different filter orders. This introduces an artificial dependence of DC group delay on the filter order. Consequently, it is important to check what type of frequency normalization is used by the chosen filter design routine.