The terms of the sequence A276123, defined by $a_0=a_1=a_2=1$ and $$a_n=\dfrac{(a_{n-1}+1)(a_{n-2}+1)}{a_{n-3}}\;,$$ are all integers (it's easy to prove that for all $n\geq2$, $a_n=\frac{9-3(-1)^n}{2}a_{n-1}-a_{n-2}-1$).
But is it also true for the sequence A276175 defined by $a_0=a_1=a_2=a_3=1$ and $$a_n=\dfrac{(a_{n-1}+1)(a_{n-2}+1)(a_{n-3}+1)}{a_{n-4}} \;\;?$$
Update : I crossposted to MO.
Yes, $(a_n)$ is a sequence of integers.
To prove this we first need to study some auxiliary sequences that satisfy a polynomial recurrence relation (unlike $(a_n)$ which has a rational fraction as its recurrence).
Consider the sequences $(b_n)$ of positive reals satisfying the recurrence relation $b_nb_{n+4} = b_{n+1}b_{n+2}b_{n+3} + 1$.
It turns out we can express $b_{n+8}$ as a polynomial in $b_n, \ldots, b_{n+7}$ :
Since $b_{n+1}b_{n+5} \equiv b_{n+2}b_{n+6} \equiv b_{n+3}b_{n+7} \equiv 1 \pmod {b_{n+4}}$ and $b_{n+1}b_{n+2}b_{n+3} \equiv -1 \pmod {b_{n+4}}$, we have $b_{n+5}b_{n+6}b_{n+7} \equiv -1 \pmod {b_{n+4}}$, which suggests the existence of a formula for $b_{n+8}$.
With this roadmap, we can write
$(b_{n+1}b_{n+2}b_{n+3})(b_{n+5}b_{n+6}b_{n+7}+1) \\ = (b_{n+1}b_{n+5})(b_{n+2}b_{n+6})(b_{n+3}b_{n+7}) + (b_{n+1}b_{n+2}b_{n+3}) \\ = (b_{n+2}b_{n+3}b_{n+4}+1)(b_{n+3}b_{n+4}b_{n+5}+1)(b_{n+4}b_{n+5}b_{n+6}+1)+(b_nb_{n+4}-1) \\ = b_{n+4}.F(b_{n+i})$
where $F$ is some big polynomial. And finally,
$(b_{n+5}b_{n+6}b_{n+7}+1) = (b_{n+5}b_{n+6}b_{n+7}+1)(b_nb_{n+4} - b_{n+1}b_{n+2}b_{n+3}) \\ = b_{n+4}(b_nb_{n+5}b_{n+6}b_{n+7}+b_n - F(b_{n+i})) = b_{n+4} G(b_{n+i})$.
And so, $b_{n+8} = G(b_{n+i})$. This means that if $b_0, \ldots, b_7 \in R$ for some subring $R$ of $\Bbb R$, then the whole sequence is in $R$.
Now to link back to the original sequence.
Given such a sequence $(b_n)$, we define a sequence $(a_n)$ by $a_n = b_nb_{n+1}b_{n+2}$.
This sequence satisfies $a_na_{n+4} = (b_n b_{n+1}b_{n+2})(b_{n+4}b_{n+5}b_{n+6}) \\ = (b_n b_{n+4})(b_{n+1} b_{n+5})(b_{n+2} b_{n+6}) = (b_{n+1}b_{n+2}b_{n+3}+1)(b_{n+2}b_{n+3}b_{n+4}+1)(b_{n+3}b_{n+4}b_{n+5}+1) \\ = (a_{n+1}+1)(a_{n+2}+1)(a_{n+3}+1)$.
Finally, taking $b_0 \ldots b_7 = \frac 12, 4, \frac 12, \frac 12, 4, \frac 12, 4, 18$, we obtain a sequence $(b_n)$ with terms in $\Bbb Z[\frac 12]$, with the corresponding $(a_n)$ sequence $1,1,1,1,8,36, \ldots$
Since the recurrence relation is symmetric, it can go backwards as well as forward, hence the ring $R_n = \Bbb Z[b_n, \ldots, b_{n+7}]$ is independant of $n$. There is no hope of finding $8$ consecutive integer values in our sequence $b_n$.
This excludes the possibility of any denominator occuring that is not a power of $2$. We will use $2$-adic analysis to exclude this last possibility below
Note that using this definition,
$a_na_{n+2}/a_{n+1}(a_{n+1}+1) = b_nb_{n+2}b_{n+4}/(b_{n+1}b_{n+2}b_{n+3}+1) = b_{n+2}$,
and so to go in the other direction you have to define $(b_n)$ from $(a_n)$ with $b_n = a_{n-2}a_n/a_{n-1}(a_{n-1}+1)$. Then, once again the recurrence relation of $(b_n)$ follows from that of $(a_n)$.
This shows that for any such rational sequence $(a_n)$, there is a corresponding rational sequence $(b_n)$, and so $(a_n)$ is in a finitely generated subring of $\Bbb Q$.
Addendum to clarify how to prove that the sequence $(a_n)$ stays in the 2-adic integers (with no deeper algebraic investigation of the sequences)
Experimenting with the recursion reveals that the behaviour of $(a_n)$ and $(b_n)$ has a pattern that approximately repeats every 11 steps. If you try to compute $a_n$ in the $2$-adics naively, no matter which transition map you use (the recursion on the $(a_n)$, on the $(b_n)$ or any polynomial map on the $(b_n)$) you end up losing precision at a steady rate.
To avoid losing precision, we have to find 2-adics neighbourhoods of the sequence where the transition maps are homeomorphisms, and even better, some iterate that are local isometries for one appropriate basis of $\Bbb Q_2^4$. To do so we have to compute the derivatives of the iterates of the transition maps and check that its characteristic polynomial has coefficient in $\Bbb Z_2$, and finally that the determinant is invertible. I will choose to work directly on the initial recurrence on $(a_n)$.
Let $f : \Bbb Q_2^4 \mapsto \Bbb Q_2^4$ be the partial map $(a_0,a_1,a_2,a_3) \mapsto (a_1,a_2,a_3,(a_1+1)(a_2+1)(a_3+1)/a_0)$.
Given a point $x = (a_0,a_1,a_2,a_3)$, assuming $f^{\circ 11}(x)$ is well-defined, we can compute its Jacobian matrix $J^{11}_x = (df_i/da_j)$, or $J^{11}_x = (da_{i+11}/da_j)$ if we consider $a_4,\ldots,a_{14}$ as functions of $a_0,a_1,a_2,a_3$.
If $\det(J^{11}_x) \neq 0$ then the iterate $f^{\circ 11}$ induces an isomorphism on the tangent spaces, and if the coefficients of the characteristic polynomial are integers, then there is a basis of the tangent space (forgetting about their different base point for a moment) where the matrix has all its entries in $\Bbb Z_2$. That is, $J^{11}_x$ is similar to a matrix in $GL_4(\Bbb Z_2)$.
As we can expect from experimentation, this does turn out to be true on points $(a_{11k},\ldots,a_{11k+3})$ of the sequence, and as luck would have it, there is a simple such basis that works simultaneously for all the points of the form $x = (a_{11k},a_{11k+1},a_{11k+2},a_{11k+3})$.
Concretely, we see the $(a_i)$ as constants, we use $f^{\circ 11}$ to map $(a_0+x_0, a_1+2x_1, a_2+2x_2, a_3+x_3)$ to $(a_{11}+y_0, a_{12}+2y_1, a_{13}+2y_2, a_{14}+y_3)$. This induces a partial map $g : \Bbb Q_2^4 \to \Bbb Q_2^4$, $(x_i) \mapsto (y_i)$ sending $0$ to $0$ and where the Jacobian at the origin $J(g)_0 \in GL_4(\Bbb Z_2)$ is an invertible matrix with $2$-adic integer coefficients.
The next step is to determine a neighbourhood $N(e) = (2^e \Bbb Z_2)^4$ of the origin where the $g|_{N(e)} : N(e) \to N(e)$ is an isometry. Here we can afford to only keep track of the order of all the higher-order terms naively. Any increase in $e$ (the lower bound on the $2$-valuation of the $x_i$) results in at least a two-fold increase in the higher-order terms, so if we choose $e$ large enough, the higher-order terms appearing in the $y_i$ will eventually stay larger then $e$. Once this happens, we have a local isometry.
To perform the operations, we use $x_i = O(2^e)$, write $\bar x$ for the $2$-valuation of $x$, and naively compute additions, multiplications, and inverses like so :
$(a + \sum a_i x_i + O(2^u)) + (b + \sum b_i x_i + O(2^v)) = ((a+b) + \sum (a_i + b_i) x_i + O(2^{\min\{u,v\}}))$
$(a + \sum a_i x_i + O(2^u)) (b + \sum b_i x_i + O(2^v)) = ((ab) + \sum (a_i b + a b_i) x_i + O(2^{\min\{ \bar{a_i}+\bar{b_j}+2e, \bar a + u, \bar b + v, \bar {a_i}+v+e, \bar {b_i}+u+e, u+v\}}))$
$(a + \sum a_i x_i + O(2^u))^{-1} = (a^{-1} - \sum a_i a^{-2} x_i + O(2^{\min \{u-2\bar a,2u-3\bar a,2\bar{a_i}+2e-3\bar a\}}))$
Using this technique we can prove that $g$ is an isometry when restricted to $N(8)$.
Then it is a matter of checking that $f^{\circ 11}$ keeps sending neighbourhoods to neighbourhoods isometrically until we get back to the starting point, which happens after $128$ steps.
After doing this and checking that the sequence stays an integer during all of the 1408 individual steps, we have proved that any sequence starting close enough to $(1,1,1,1)$ will generate a sequence of $2$-adic integers.
The size $2^{8}$ for my neighbourhoods may not be optimal, and the amount of computation needed to get back to the starting point only depends linearly on the size, so it's not that important to get the best possible.
The observation that we repeat so quickly suggests that there are 3 hidden invariants who cut out the curves that each sequence traces inside the neighbourhoods : $a_{11k}$ modulo $2^n$ determines $a_{11k+1}$ and $a_{11k+2}$ modulo $2^{n+1}$ and $a_{11k+3}$ modulo $2^n$
Additionally, when I looked at polynomial maps on octuplets $(b_0,\ldots,b_7) \mapsto (b_{11},\ldots,b_{18})$ instead, I found that the Jacobians were $x^8+x^4 = (x+1)^4x^4$ modulo $2$. I didn't look for a nice basis there, but having 4 additional even eigenvalues should imply that if you start with a nearby random octuplet, it should get attracted closer and closer to the octuplet of a "real" sequence of b_n generated by 4 inputs and the original relation. So those polynomial maps are also remarkably stable.
Following up on a comment by David E. Speyer (indicating that the $b_n$ sequence arises from the cluster algebra $\mathcal{A}(Q)$ associated to a certain quiver $Q$), I would like to mention that the $a_n$ sequence arises from the Poisson cluster variety $\mathcal{X}(Q)$ associated to the same quiver $Q$. Moreover, the relation $$a_n = b_n b_{n+1} b_{n+2}$$ explained in mercio's answer is an instance of (the pullback of) Fock-Goncharov's $\textit{cluster ensemble}$ map $p: \mathcal{A}(Q) \to \mathcal{X}(Q)$. This point of view can be applied to show that, for instance, the sequence $(c_n)$ defined by $c_0 = 1, c_1 = 2, c_2 = 6, c_3 = 42$ and for $n \geq 4$, $$c_n = \frac{(1+c_{n-1})(1+c_{n-2})(1+c_{n-3})}{c_{n-4}},$$ consists entirely of (positive) integers (the $c_n$ sequence consists of values in a so-called Y-frieze pattern associated to $Q$).
Perhaps (with a bit more work) this point of view could help settle (in the affirmative) the OP's question.