Does $\gcd$ of values of $n$-th degree primitive integer polynomial always divide $n!$? Let's consider the following conjecture:

Let $f(x)=a_0+a_1x+a_2x^2+\dots +a_nx^n$ be a primitive non-zero polynomial with integer coefficients and let $g=\gcd(\{f(k), k \in \mathbb{Z}\})$, then $g$ divides $n!$.

Or equivalently just:

Let $0\not\equiv f(x)\in \mathbb{Z}[x]$, then $\gcd(f(k)) \mid c(f)\cdot (\deg f)!$.

Here $c(f)$ is the content of a polynomial.
Can we prove/disprove this?
Edit: Found literature reference
Today I found that $g$ is called the fixed divisor of a polynomial and with this I was able to trace the same claim to the article Uber ganzwertige ganze funktionen by George Pólya from 1915. Unfortunately the article is in german so I am not sure how the proof goes (various articles mention this claim and cite this one without providing the proof themselves).
Some thoughts
I found this based on observation for small degrees. For $n=2$ and $\gcd(a_0,a_1,a_2)=1$ we can prove $g \mid 2$ (see below), and similarly for $n=3$ we can show $g \mid 6=3!$. I also wrote a simple script in Maple to verify over small degrees and coefficients range and haven't found a counterexample yet, so I wonder if this can be proven to hold in general.
It can be useful that $\gcd(\{f(k), k \in \mathbb{Z}\})=\gcd(f(0),f(1),\dots,f(n))$ (or any consecutive $n+1$ integers for that matter). Currently I am thinking that Lagrange interpolation could be perhaps useful in generic case, I worked out that if we know $f$ on consecutive integers $i=0,1,\dots,n$, then the Lagrange interpolation implies
$$
f(x)=\sum_{i=0}^{n}  (-1)^{n-i}\binom{x}{i}\binom{x-i-1}{n-i}f(i).
$$
So I guess if we assume $\gcd(f(i))=d$ where $d \not\mid n!$, we should be able to use the above to reach contradiction with $f(x)$ having integer and coprime coefficients, but I wasn't able to.
Proof for $n=2$
For $f(x)=a_0+a_1x+a_2x^2$ and $\gcd(a_0,a_1,a_2)=1$, we want to show $\gcd(f(-1),f(0),f(1)) \mid 2$.
Now $f(-1)=a_0-a_1+a_2$, $f(0)=a_0$, $f(1)=a_0+a_1+a_2$ and so
\begin{align}
\gcd(f(-1),f(0),f(1))
&=\gcd(a_0-a_1+a_2,a_0+a_1+a_2,a_0)\\
&=\gcd(-a_1+a_2,a_1+a_2,a_0)\\
&=\gcd(2a_2,a_1+a_2,a_0)\\
\end{align}
Now $g=\gcd(2a_2,a_1+a_2,a_0) \mid 2a_2$, but $a_2$ is coprime to $g$ (otherwise we had a prime $p\mid a_2, p\mid a_1+a_2, p \mid a_0$ and so $p \mid \gcd(a_2,a_1,a_0)=1$, impossible). Hence $\gcd(2a_2,a_1+a_2,a_0) \mid 2$.
 A: Here is a proof using $g:=\gcd\{f(d), d \geq 1\}=\gcd(b_0,b_1,\dots,b_n)$ where $f(x)$ is expressed as
$$f(x)=b_0+b_1\binom{x}{1}+\dots+b_n\binom{x}{n}.$$
(This is a known result, see this my older post for example)
Furthermore, having $f(x)=a_0+a_1x+a_2x^2+\dots +a_nx^n$, we can express $b_i$'s by
$$
b_i = i!\sum_{k=i}^{n}{k\brace i}a_k.
$$
Proof: This follows directly from a well known identity for Stirling numbers of the second kind:
$$\sum_{i=0}^{k}{k\brace i}(x)_i=x^k.$$
We have $(x)_i=x(x-1)\cdots(x-i+1)=\binom{x}{i}i!$ and then $$f(x)=\sum_{k=0}^{n} a_k x^k=\sum_{k=0}^{n} a_k \sum_{i=0}^{k}{k\brace i}\binom{x}{i}i!=\sum_{i=0}^{n}\binom{x}{i}\sum_{k=i}^{n} a_k {k\brace i}i!.$$
$\square$
Hence $g=\gcd(b_0,b_1,\dots,b_n)$ and so $g \mid b_i = i!\sum_{k=i}^{n}{k\brace i}a_k$. Now we prove $g\mid n!a_i$ for all $i=0,1,\dots,n$ by induction. First notice that $g\mid b_n= n!\sum_{k=n}^n {k\brace n}a_k=n!a_n$ (a base case). Now assume $g\mid n!a_n, g\mid n!a_{n-1},\dots g\mid n!a_{n-i}$, then
\begin{align}
g 
&\mid b_{n-i-1}=(n-i-1)!\sum_{k=n-i-1}^{n}{k\brace n-i-1}a_k\\
\\&\mid n!\sum_{k=n-i-1}^{n}{k\brace n-i-1}a_k\\
\\&= n!{n-i-1\brace n-i-1}a_{n-i-1}+ n!\sum_{k=n-i}^{n}{k\brace n-i-1}a_k\\
\\&= n!a_{n-i-1}+ \sum_{k=n-i}^{n}{k\brace n-i-1}n!a_k\\
\end{align}
and since by the induction hypothesis $g\mid n!a_k$ in the sum, it follows $g \mid n!a_{n-i-1}$.
So finally we have $g\mid n!a_i$ for $i=0,1,\dots,n$ and hence
$$
g\mid (n!a_0,n!a_1,\dots,n!a_n)=n!\gcd(a_0,a_1,\dots,a_n).
$$
Particularly for a primitive polynomial $\gcd(a_0,a_1,\dots,a_n)=1$ and
$$
g \mid n!.
$$
A: The equations $f(k)=a_0+a_1k+\dots+a_nk^n$ for $k=1,\dots,n+1$ give us $f(k)$'s as an integer linear combination of $a_0,a_1,\dots,a_n$. In a matrix notation we have $V_n=(i^{j-1})$ and
$$
\begin{pmatrix}
f(1)\\
f(2)\\
\vdots\\
f(n+1)
\end{pmatrix}=V_n\begin{pmatrix}
a_0\\
a_1\\
\vdots\\
a_n
\end{pmatrix}.
$$
However per this post and especially the method referred by Jean Marie it follows that
$$
n!\begin{pmatrix}
a_0\\
a_1\\
\vdots\\
a_n
\end{pmatrix}=n!V_n^{-1}\begin{pmatrix}
f(1)\\
f(2)\\
\vdots\\
f(n+1)
\end{pmatrix}
$$
where elements of matrix $n!V_{n}^{-1}$ are
$$
(-1)^{i+j}\sum_{m=i}^{n+1}\frac{n!}{(m-1)!}{m\brack i}\binom{m-1}{j-1}
$$
which are clearly integers. Hence
$$
n!a_i=b_{i,1}f(1)+b_{i,2}f(2)+\dots+b_{i,n}f(n)
$$
for some integers $b_{i,j}$, and since $g\mid f(k)$, we have $g\mid n!a_i$ and so
$$
g\mid \gcd(n!a_0,n!a_1,\dots,n!a_n)=n!\gcd(a_0,a_1,\dots,a_n).
$$
