Thank you for posting this question, I enjoyed trying to answer it.
Start with the expression that Mathematica gave you and replace each
argument $\frac14$ of a hypergeometric function with $\frac z4$,
because we will be taking limits. I will call the two hypergeometric
functions $Q_1(z)$ and $Q_2(z)$. Each term can be brought to a closed
form by using identity 16.6.2 from the DLMF.
Setting $a=\frac12$, $b=1-\frac\nu2$, we get
$$ Q_1(z) = (1-z)^{-\frac12} F\left(\frac16,
\frac36,\frac 56; 1-\frac\nu2, 1+\frac\nu2; \frac{-27 z}{4(1-z)^3}
\right), $$
and setting $a=\frac{1+3\nu}{2}$, $b=1+\frac\nu2$, we get
$$ Q_2(z) = (1-z)^{-\frac{1+3\nu}{2}} F\left(\frac{1+3\nu}6,
\frac{3+3\nu}{6}, \frac{5+3\nu}{6}; 1+\frac\nu2, 1+\nu;
\frac{-27z}{4(1-z)^3} \right). $$
(Note that there are 6 possible identities to try per function, one for each possible choice of
$a$ and $b$ from the parameters, so it helps to do this on a
computer.)
The reason this works is that now the point $z=1$ is a singular point
of the hypergeometric functions on the right hand side, and Mathematica
will succeed in finding the limits as $z\to1$. The expression for the
whole integral that you have is
$$ Q =
\frac{2^{\frac43}\pi^{\frac12}}{\Gamma(-\frac16)\Gamma(\frac56+\frac\nu2)\Gamma(\frac56-\frac\nu2)\sin\frac{\nu\pi}2}
\left( -1 + 3^{-\frac{3\nu}2}\cos\left(\frac{\nu\pi}{2}\right)
\frac{\Gamma(\frac{1+3\nu}{2})
\Gamma(\frac56-\frac\nu2)}{\Gamma(\frac{1+\nu}2)\Gamma(\frac56+\frac\nu2)}
\right). $$
Call the large expression in brackets $A$, and then write
$$ A = -1 + B 3^{-\frac{3\nu}{2}}\cos\frac{\pi\nu}{2}, \qquad B = \frac{\Gamma(\frac{1+3\nu}{2})
\Gamma(\frac56-\frac\nu2)}{\Gamma(\frac{1+\nu}2)\Gamma(\frac56+\frac\nu2)}. $$
Now, Mathematica will not simplify $A$ or $B$ on its own, so it needs
help. Set $x=\frac16+\frac\nu2$, and use the multiplication formula to
get
$$
\frac{\Gamma(\frac{1+3\nu}2)}{\Gamma(\frac{1+\nu}{2})\Gamma(\frac56+\frac\nu2)}
= \frac{\Gamma(3x)}{\Gamma(x+\frac13)\Gamma(x+\frac23)} =
\frac{\Gamma(x)}{2\pi} 3^{3x-1/2}. $$
After this, $A$ simplifies to
$$ A = -1 +
\frac{\Gamma(\frac16+\frac\nu2)\Gamma(\frac56-\frac\nu2)}{2\pi}\cos\frac{\pi\nu}{2}
= -1 + \frac{\cos\frac{\pi\nu}{2}}{2\sin(\frac\pi6+\frac{\pi\nu}{2})},
$$
where I've also used the reflection formula for $\Gamma(z)\Gamma(1-z)$ to get rid of the
gamma functions. Some further amount of manual trigonometry yields
$$ A = -\frac{\sqrt{3}}{2}\frac{\sin\frac{\pi\nu}{2}}{\sin(\frac\pi6 +
\frac{\nu\pi}{2})}. $$
Finally, write
$$ \frac{1}{\sin(\frac\pi6+\frac{\pi\nu}{2})} =
\frac{\Gamma(\frac16+\frac\nu2)\Gamma(\frac56-\frac\nu2)}{\pi}, $$
and substitute back. Lots of things cancel, and the answer is
$$ Q = -\frac{3^{1/2}2^{1/3}}{\pi^{1/2}}
\frac{\Gamma(\frac16+\frac\nu2)}{\Gamma(-\frac16)\Gamma(\frac56+\frac\nu2)}. $$
This closed form is equivalent to the one you gave through the
use of $\Gamma(\frac16)\Gamma(-\frac16)=-12\pi$.
P.S.
I would also like to note that the integral
$$ I(\nu,c) = \int_0^\infty J_\nu(x)^2 J_\nu(c x)\,dx $$
and its general form
$$ \int_0^\infty x^{\rho-1}J_\nu(a x) J_\mu(b x) J_\lambda(c x)\,dx $$
appear in Gradshteyn and Ryzhik, and you can find a paper "Some
infinite integrals involving bessel functions, I and II" by
W. N. Bailey, which evaluates this integral in terms of Appell
functions, but only in the case $c>2$ ($|c|>|a|+|b|$), which is where the $F_4$ Appell
function converges. DLMF 16.16.6 actually gives a way to write this
integral as
$$
\frac{\Gamma(\frac{1+3\nu}{2})c^{-1-2\nu}}{\Gamma(1+\nu)^2\Gamma(\frac{1-\nu}{2})}
\,\,\,{}_2F_1\left( \frac{1+\nu}{2}, \frac{1+3\nu}{2}; 1+\nu; x
\right)^2, \qquad x = \frac{1-\sqrt{1-4/c^2}}{2}, $$
but the issue is that this is only correct for $c>2$, and the rhs is
complex for $c<2$. Appell function would only be defined by analytic
continuation in this case anyway, and I didn't find anything useful about
non-principal branches of Appell or hypergeometric functions.
For $c>2$, Mathematica also gives the following:
$$
I(\nu,c) = \frac{\Gamma(\frac{1+3\nu}{2})c^{-1-2\nu}}{\Gamma(\frac{1-\nu}{2})\Gamma(1+\nu)^2}
\,\,\,{}_3F_2\left( \frac{1+\nu}{2}, \frac{1}{2}+\nu,
\frac{1+3\nu}{2}; 1+\nu, 1+2\nu; \frac{4}{c^2} \right), $$
but this is incorrect when $c<2$.
For the first identity, you can show that your expression for the beta function is equivalent to the integral formula
$$ B(x,y) = \int_0^1 t^{x-1} (1-t)^{y-1} dt $$
Then expand the second factor in the integrand,
$$ B(x,y) = \int_0^1 t^{x-1} \sum_{k=0}^{\infty} (-1)^k \left(
\begin{array}{c}
y-1 \\
k
\end{array}
\right) t^k dt $$
Given that each term of the summand is measurable we can use apply the monotone convergence theorem and switch the orders of integration and summation to give
$$ B(x,y) = \sum_{k=0}^{\infty} (-1)^k \left(
\begin{array}{c}
y-1 \\
k
\end{array}
\right) \int_0^1 t^{x+k-1} dt $$
$$ = \sum_{k=0}^\infty \left(
\begin{array}{c}
k-y \\
k
\end{array}
\right) \frac{1}{x+n} $$
where in the last line we used an identity for the binomial coefficient and computed the integral.
The second result is merely a restatement of the first, writing the binomial coefficient in terms of Gamma functions.
Best Answer
This is not a complete answer. Using the limit representation in Tyma Gaidash's answer, and the known results $$ H_{n + j + r} = \log (n + j + r) + \gamma + o(1), \quad r\to+\infty, $$ (with $\gamma$ being the Euler–Mascheroni constant) and $$ \sum\limits_{j = 0}^n {\binom{n}{j}^2 } = \binom{2n}{n}, $$ we can derive the explicit formula $$ \boxed{I_n = \binom{2n}{n}\gamma - \sum\limits_{j = 0}^n {\binom{n}{j}^2 (2(H_{n - j} - H_j )\log (( n+j)!) + H_{n + j} )} ,} $$ provided that $$ \mathop {\lim }\limits_{r \to + \infty } \sum\limits_{j = 0}^n {\binom{n}{j}^2 \left( {2(H_{n - j} - H_j )\log ((n + j + r)!) + \log (n + j + r)} \right)} = 0. $$ Denote the expression under the limit by $S_n(r)$. Performing the change of summation index $j\to n-j$ and taking the average with the original expression, we find $$ S_n (r) = \sum\limits_{j = 0}^n {\binom{n}{j}^2 \left( {(H_{n - j} - H_j )\log \frac{{(n + j + r)!}}{{(2n - j + r)!}} + \frac{1}{2}\log ((2n - j + r)(n + j + r))} \right)} . $$ Stirling's formula and the Maclaurin series of the logarithm then yields $$ S_n (r) = (\log r)\sum\limits_{j = 0}^n {\binom{n}{j}^2 \left( {(H_{n - j} - H_j )(2j - n) + 1} \right)} + o(1) $$ as $r\to+\infty$. Therefore, it remains to show that $$ \sum\limits_{j = 0}^n {\binom{n}{j}^2 ((H_{n - j} - H_j )(2j - n) + 1)} =0 $$ for all $n\ge 1$. Because of the symmetry in $H_{n - j} - H_j$, this may be further simplifed to the claim that $$ \sum\limits_{j = 0}^n {\binom{n}{j}^2 (2j(H_{n - j} - H_j )+ 1)} =0 $$ for all $n\ge 1$. Computer algebra software confirms this for $n=1,2,3\ldots,200$.
Addendum. The limit expression in the original version of Tyma Gaidash's answer was $$ I_n=\lim_{r\to\infty}\sum_{j=0}^n\binom n j^2\left(2 (H_{n-j}-H_j)\ln\left(\frac{(n+j+r)!}{(n+j)!}\right)+H_{n+j+r}-H_{n+j}\right). $$