How can it be proved that the sequence $(a_{n})_{n\geq0}$ satisfying the recurrence $$(n+1)a_{n} - r(r-1)(r(n-1)+1)a_{n-1}=0\quad a_{0}=1 \quad r \geq 2 \in \mathbb{N}$$ is always a sequence of integers for any $r$?
In this case, the recurrence for the Catalan numbers (which are integers) is given by $r=2$: $(n+1)a_{n} - 2(2(n-1)+1)a_{n-1}=(n+1)a_{n} - (4n-2)a_{n-1}=0\quad a_{0}=1$
As it was pointed out in the commentary $$a_{n}=\dfrac{((r - 1) r^2)^n Γ\left(n + \dfrac{1}{r}\right)}{Γ(n + 2) Γ\left(\dfrac{1}{r}\right)}$$ Using the formula found in https://en.wikipedia.org/wiki/Particular_values_of_the_gamma_function $$Γ\left(n + \dfrac{1}{r}\right)=Γ\left(\dfrac{1}{r}\right)\dfrac{(rn-r+1)!^{(r)}}{r^n}$$ where $(n)!^{(r)}$ is the $r^{th}$ multifactorial of $n$.
Therefore, for $n\in \mathbb{N}_{0}$ and $r\geq 2 \in \mathbb{N}$
$$a_{n}=\dfrac{(r-1)^{n}r^{n}(r(n-1)+1)!^{(r)}}{(n+1)!}$$
Thus, the problem reduces to show that $(n+1)!$ divides $(r-1)^{n}r^{n}(r(n-1)+1)!^{(r)}$
The motivation behind the question arises from a computational experiment to find rational approximations $\dfrac{p_{n}}{q_{n}}$ to $Γ\left(\dfrac{1}{r}\right)$, looking for a proof of the irrationality of this number.
Update: The following formula can be generalized to higher $r$, but we will deal with $r=3$ for simplicity:
$$a_{n}=\dfrac{2^{n}3^{n}(3n-2)!^{(3)}}{(n+1)!}$$
By the definition of the multifactorial:
$$(3n-2)!^{(3)}(3n-1)!^{(3)}(3n)!^{(3)}=(3n)!$$
and
$$(3n)!^{(3)}=\prod_{i=1}^{n} 3i=3^{n}n!$$
So,
$$3^n(3n-2)!^{(3)}=\dfrac{1}{(3n-1)!^{(3)}}\dfrac{3n!}{n!}$$
Plugging this into the first expression
$$a_{n}=\dfrac{2^{n}}{(3n-1)!^{(3)}}\dfrac{3n!}{n!}\dfrac{1}{(n+1)!}$$
and multiplying many times by 1 ($2n!/2n!$, $n!/n!$), we get
$$a_{n}=\dfrac{2^{n}n!}{(3n-1)!^{(3)}}\dfrac{3n!}{n!2n!}\dfrac{2n!}{(n+1)!n!}=\dfrac{2^{n}n!}{(3n-1)!^{(3)}}\binom{3n}{n}\dfrac{1}{n+1}\binom{2n}{n}=\dfrac{2^{n}n!}{(3n-1)!^{(3)}}\binom{3n}{n}C_{n}$$
where $C_{n}$ are the Catalan numbers
So now the main question is: what is the $p$-adic valuation (exponent of the prime $p$ in the factorization of the number) of the expression on the left?
$$\nu_{p}\left(\dfrac{2^{n}n!}{(3n-1)!^{(3)}}\right)$$
The primes whose $p$-adic valuation is negative will be in the denominator, and these should divide $\binom{3n}{n}C_{n}$. If we know this, in theory, we should be able to use Kummer's theorem to settle the divisibility question, but I wasn't able to tackle it.
Interestingly, while trying some values of $n$, not all of these prime factors that are in the denominator divide $\binom{3n}{n}$ : some of them divides $\binom{3n}{n}$ and some of them divides $C_{n}$.
Example with $n=13$:
$$\dfrac{2^{13}13!}{(3*13-1)!^{(3)}}\binom{3*13}{13}=\dfrac{15412156416}{115}$$
but $115=5*23|C_{13}$
We first show that, for every prime $p\nmid r$, $$\nu_p((n+1)!)\leq \sum_{i=-1}^{n-1}\nu_p(ri+1),$$ where $\nu_p$ for negative integers is defined exactly the same as for positive integers. By writing $$\nu_p(x)=\sum_{k=0}^\infty [1\text{ if }p^k|x],$$ it suffices to show that $$\left\lfloor\frac{n+1}{p^k}\right\rfloor$$ is at most the number of $-1\leq i<n$ for which $p^k|ri+1$. This is equivalent to the number of $-1\leq i<n$ lying in a particular residue class modulo $p^k$. In a sequence of $p^k$ consecutive integers, as $\gcd(p^k,r)=1$, there will be exactly one integer $i$ with $p^k|i$; as a result, there are at least $\lfloor (n+1)/p^k\rfloor$ integers $i$ in the range $[-1,n)$ (of length $n+1$) that are multiples of $p^k$.
In addition, for every prime $p|r$, $\nu_p((n+1)!)\leq n\leq \nu_p(r^n)$, with equality only liable to hold when $p=2$ and $r\equiv 2\pmod 4$.
So, $$\frac{r^n(1-r)(1)(1+r)(1+2r)\cdots(1+(n-1)r)}{(n+1)!}$$ is an integer. Your $a_n$ is $-(r-1)^{n-1}$ times this quantity, so it is an integer as well.