On $\int_0^1\arctan\,_6F_5\left(\frac17,\frac27,\frac37,\frac47,\frac57,\frac67;\,\frac26,\frac36,\frac46,\frac56,\frac76;\frac{n}{6^6}\,x\right)\,dx$

I think you should start by the following general form

$${}_{k}F_{k-1}\left(\frac{1}{k+1} ,\cdots ,\frac{k}{k+1};\frac{2}{k} \cdots ,\frac{k-1}{k},\frac{k+1}{k};\left( \frac{m(1-m^k)}{f_k}\right)^k \right) = \frac{1}{1-m^k}$$

Where

$$f_k \equiv \frac{k}{(1+k)^{1+1/k}}$$

Put $k=4$

$${}_{4}F_{3}\left(\frac{1}{5} ,\cdots ,\frac{4}{5};\frac{2}{4} \cdots ,\frac{3}{4},\frac{5}{4};\left( \frac{m(1-m^4)}{f_4}\right)^4 \right) = \frac{1}{1-m^4}$$

The argument simplifies to

$$ \left(\frac{m(1-m^4)}{\frac{4}{5^{1+1/4}}}\right)^4 = \frac{5^5 \,m^4(1-m^4)^4}{4^4}$$

Hence we have

$${}_{4}F_{3}\left(\frac{1}{5} ,\cdots ,\frac{4}{5};\frac{2}{4} \cdots ,\frac{3}{4},\frac{5}{4};\frac{5^5 (1-m)m^4}{4^4} \right) = \frac{1}{m}$$

Now suppose we want to find

$$\int \arctan\,_4F_3\left(\frac15,\frac25,\frac35,\frac45;\,\frac24,\frac34,\frac54;\,\frac{n}{4^4}\,x\right)\,\mathrm dx$$

Use $nx = 5^5m^4(1-m)$

$$\frac{5^5}{n}\int (4m^3(1-m)-m^ 4)\arctan\,_4F_3\left(\frac15,\frac25,\frac35,\frac45;\,\frac24,\frac34,\frac54;\,\frac{5^5 (1-m)m^4}{4^4}\right)\,\mathrm dm$$

This simplifies to

$$\frac{5^5}{n}\int (4m^3(1-m)-m^ 4)\arctan\,\frac{1}{m}\,\mathrm dm$$

The anti-derivative of this is elementary.

Note that when we use substittution we will need the roots of

$$m^5-m^4+\frac{n}{5^5}= 0$$

As conjectured by the OP.

Using the value $n=4$ we can verify Reshetnikov result

\begin{align}\int_0^1\arctan{_4F_3}\left(\frac15,\frac25,\frac35,\frac45;\frac24,\frac34,\frac54;\frac{1}{64}\,x\right)\,dx &= \frac{5^5}{4}\int^{\alpha}_1 (4x^3(1-x)-x^ 4)\arctan \left(\frac{1}{x} \right)\,\mathrm dx\\&=0.7857194\dots \end{align}

Where $\alpha$ is the largest positive root of

$$m^5-m^4+\frac{4}{5^5}= 0$$


The result does generalize. We will derive a closed formula for the integral $$I_k(n) = \int_0^1\arctan\left[_{k}F_{k-1}\left(\frac{1}{k+1},\cdots, \frac{k}{k+1},\frac{2}{k},\cdots,\frac{k+1}{k}; \frac{nx}{k^k}\right)\right]\,{\rm d}x$$ valid for $|n| \leq k^k$.

For this particular combination of arguments the hypergeometrical function simplifies greatly, see property $(25)$ here$^{(1)}$. Using this we can perform a substitution $\frac{nx}{(k+1)^{k+1}} = y^k(1-y^k)^k$ and the integrand simplifies to $\frac{(k+1)^{k+1}}{n}\frac{1}{1-y^k}\frac{d}{dy}\left(y^k(1-y^k)^k\right)$. Now performing integration by parts and the substitutions $y^k\to y \to 1-y$ we arrive at

$$I_k(n) = \frac{\pi}{2} - \arctan(z_0) + \frac{(k+1)^{k+1}}{n}\int_{z_0}^1\frac{w^{k+1}-w^k}{1+w^2}{\rm d}w$$

where $z_0$ is the largest positive root of $z_0^{k+1} - z_0^k + \frac{n}{(k+1)^{k+1}} = 0$. The integral above can be evaluated in closed form using

$$\matrix{\int \frac{w^{2k}}{1+w^2}{\rm d}w &=& (-1)^k\left[\arctan(w) + F_k(w)\right] + C\\\int \frac{w^{2k+1}}{1+w^2}{\rm d}w &=& (-1)^k\left[\frac{1}{2}\log(1+w^2) + G_k(w)\right] + C}$$ where $F_k(x) = \sum_{i=1}^{k}\frac{(-1)^i x^{2i-1}}{2i-1}$ and $G_k(x) = \sum_{i=1}^{k}\frac{(-1)^i x^{2i}}{2i}$. These relations are easily proven by differentation and summing a geometrical series. This leads to

$$I_{2k}(n) = \frac{\pi}{2} - \arctan(z_0) + \frac{(-1)^k(2k+1)^{2k+1}}{n}\left[\frac{1}{2}\log\left(\frac{2}{1+z_0^2}\right) - \frac{\pi}{4} + \arctan(z_0)\right] + \frac{(-1)^k(2k+1)^{2k+1}}{n}\left[G_k(1) - G_k(z_0) - F_k(1) + F_k(z_0)\right]$$ $$I_{2k+1}(n) = \frac{\pi}{2} - \arctan(z_0) + \frac{(-1)^{k}(2k+2)^{2k+2}}{n}\left[-\frac{1}{2}\log\left(\frac{2}{1+z_0^2}\right) - \frac{\pi}{4} + \arctan(z_0)\right] + \frac{(-1)^{k}(2k+2)^{2k+2}}{n}\left[- G_k(1) + G_k(z_0) -F_{k+1}(1) + F_{k+1}(z_0)\right]$$

For the particular case $k = 2$ ($p=3$ in the question) you asked for this reduces to

$$I_2(1) = \frac{27}{2}\log \left(\frac{z_0^2+1}{2}\right) - \frac{27}{2}(z_0-1)^2 - 28\arctan(z_0) + \frac{29 \pi}{4} \simeq 0.795296 \\ \text{where } z_0 = \frac{1}{3} \left(1+2 \cos \left(\frac{\pi }{9}\right)\right)$$

Finlly since the expression has many terms we should check that it gives the correct result. Here is some code used to check this in Mathematica:

(* Working precision *)
ndigits = 100;

Block[{$MinPrecision = ndigits},

  (* Pick k and n *)
  k = 4; n = 1;

  (* Solve for z0 *)
  z0 = Max[N[x /. Solve[{x^(k + 1) - x^k + n/(k + 1)^(k + 1) == 0, x > 0}, x], ndigits]];

  (* Define the integral *)
  ak = Table[i/(k + 1), {i, 1, k}];
  bk = Table[If[i == k, (i + 1)/k, i/k], {i, 2, k}];
  exact = NIntegrate[ArcTan[HypergeometricPFQ[ak, bk, (x n)/k^k]], {x, 0, 1}, WorkingPrecision -> ndigits];

  (* Our solution *)
  result = If[Mod[k, 2] == 0,

     (* For even k *)
     \[Pi]/2 - ArcTan[z0] + ((-1)^(k/2) (k + 1)^(k + 1))/ n (1/2 Log[2/(1 + z0^2)] - \[Pi]/4 + ArcTan[z0] + Sum[(-1)^i ((1 - z0^(2 i))/(2 i) - (1 - z0^(2 i - 1))/(2 i - 1)), {i, 1, k/2}])

     ,

     (* For odd k *)
     \[Pi]/2 - ArcTan[z0] + ((-1)^((k - 1)/2) (k + 1)^(k + 1))/n (-(1/2) Log[2/(1 + z0^2)] - \[Pi]/4 + ArcTan[z0] - Sum[(-1)^i ((1 - z0^(2 i))/(2 i)If[i == (k - 1)/2 + 1, 0, 1] + (1 - z0^(2 i - 1))/(2 i - 1)), {i, 1, (k - 1)/2 + 1}])

  ];

  (* Difference between integral and our formula *)
  Print["Difference: ", N[result - exact]];

];

This demonstrates that we get the same result (to numerical precision) for a wide range of $k$ and $n$.

$^{(1)}$ It's easy to check that this holds to arbitrary precisison (checked to a few hundred digits) at least for the smallest values of $k$. Marty Cohen found this reference discussing a very similar type of property which might be useful if anyone wants to try to derive it (which likely is a very messy affair).