Let $b_{(i,j)}$ for $1 \le i, j \le n+1$ be positive integers in which for every $i, j, k, l$ where $\{i, j\} \ne \{k, l\}$, we have $(b_{(i,j)}, b_{(k,l)}) = 1$ and $b_{(i,j)} = b_{(j,i)}$. According to Lagrange interpolation formula, there exists a unique polynomial $P(x)$ with rational coefficients of degree at most $n$ that satisfies $\forall 1 \le i \le n+1, P(a_i) = \prod_{j=1}^{n+1} b_{(i,j)}$. Then we'll have:
$$(P(a_i), P(a_j)) = b_{(i,j)} > 1, \quad (P(a_i), P(a_j), P(a_k)) = (b_{(i,j)}, P(a_k)) = 1$$
Therefore, it suffices to prove that there is a sequence $b_{(i,j)}$ in which the resulting polynomial $P$ has integer coefficients in lieu of rational. According to Lagrange interpolation formula,
$$P(x) = \sum_{i=1}^{n} P(a_i) \frac{\prod_{j=1, j \neq i}^{n} (x - a_j)}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)}$$
Define $c$ to be $c = \prod_{i=1}^{n+1} \prod_{j=1, j \neq i}^{n+1} (a_i - a_j)$. Redefine $b_{(i,j)}$ so that in addition to former conditions, for every $(i, j)$, $b_{(i,j)} \equiv 1 \pmod{c}$. Obviously such sequence $b_{(i,j)}$ exists, as it's equivalent to existence of $\frac{n(n+1)}{2}$ pairwise distinct positive integers all being $1$ modulo $c$, and the sequence
$$c_1 = 1, \quad \forall 2 \le i \le \frac{n(n+1)}{2} : c_i = c \prod_{j=1}^{i-1} c_j + 1$$
works well. Furthermore, $P(a_i) = \prod_{j=1}^{n+1} b_{(i,j)} \equiv 1 \pmod{c}$.
It remains to prove that for such $P(a_i)$, $P$ would have integer coefficients. Define $k_j$ to be $k_j = \frac{P(a_j)-1}{c}$. According to Lagrange interpolation formula for constant polynomial $1$, we have
$$\begin{align*} 1 &= \sum_{i=1}^{n} \frac{1 \cdot \prod_{j=1, j \neq i}^{n} (x - a_j)}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)} = \sum_{i=1}^{n} \frac{(P(a_i) - c k_i) \prod_{j=1, j \neq i}^{n} (x - a_j)}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)} \\ &= \sum_{i=1}^{n} \frac{(P(a_i)) \prod_{j=1, j \neq i}^{n} (x - a_j)}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)} - \sum_{i=1}^{n} \frac{c k_i \prod_{j=1, j \neq i}^{n} (x - a_j)}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)} \\ &= P(x) - \sum_{i=1}^{n} k_i \prod_{j=1, j \neq i}^{n} (x - a_j) \frac{c}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)} \end{align*}$$
But by definition $\prod_{j=1, j \neq i}^{n} (a_i - a_j) \mid c$. Therefore,
$$\sum_{i=1}^{n} k_i \prod_{j=1, j \neq i}^{n} (x - a_j) \frac{c}{\prod_{j=1, j \neq i}^{n} (a_i - a_j)}$$
has integer coefficients, and so does $P(x)$.