The same Puiseux-series tactic works (effectively) to list all $y$
for which $p(X,y)$ has a quadratic factor. For example, the linear
coefficient of such a factor is $a = -(r_i+r_j)$ for some distinct $i,j$
where $r_1,\ldots,r_5$ are the roots of $p(X,y)$. But for large $y$
there are three real and two imaginary roots, and thus four choices of
$\{i,j\}$ that make $r_i+r_j$ real; and in each case you can expand
$a$ in a Laurent series in $y^{-1/6}$ with rational coefficients.
You can then concoct a polynomial in $a$ and $y$ with integer coefficients
that decays as a positive power of $1/y$ but is not identically zero,
and thus cannot be an integer once $y$ is large enough. Surely your computation
up to $10^9$ will then be way more than enough to establish the desired result.
[Added later] Here are explicit power series and polynomials
that can be used for this purpose. Let $z = y^{-1/6}$. Then
for large $y>0$ (i.e. small $z>0$) the roots of $p(X,y)$ are
$$
r_1, r_2 = \pm z^{-3} - z^6 \mp z^9 \pm \frac32 z^{15} + 5 z^{18} + O(z^{21}),
$$ $$
r_3 = z^{-4} + \frac13 z^{-2} - \frac1{81} z^2 + \frac{82}{243} z^4
+ \frac23 z^6 + \frac{5585}{6561} z^8 + O(z^{10}),
$$
and complex conjugate roots $r_4,r_5$ obtained from $r_3$ by
multiplying $z$ by a cube root of unity. Let's warm up by showing
that for large $y$ there are no integral roots, i.e. none of
$r_1,r_2,r_3$ can be an integer (we needn't worry about $r_4,r_5$
because they're not even real). For $x=r_1,r_2$ it is enough to observe that
$$
x^2 - y = \mp 2 z^3 - 2 z^6 + 4 z^{12} \pm 12 z^{15} + 12 z^{18} + O(z^{21}),
$$
which for large $y$ is a small nonzero real number and thus cannot be
an integer. For $x=r_3$ we use the $11$ monomials in $x,y$
that have a pole of order at most $9$ in $z^2$, and find the polynomial
$$
y^3 + (3-2x) y^2 + (-x^3-3x^2+4x-4) y + 3x^4-3x^3-x^2-7x+4
$$
with power series
$$
4 z^4 - \frac43 z^6 - \frac{104}{9} z^8 - \frac{980}{81} z^{10}
- \frac{32}{3} z^{12} + O(z^{14}),
$$
which again is nonzero but small for large $y$.
For a quadratic factor the candidates are
$(X-r_1)(X-r_2)$, $(X-r_1)(X-r_3)$, $(X-r_2)(X-r_3)$, and $(X-r_4)(X-r_5)$.
The first of these has $r_1+r_2 = -2z^6 + O(z^{18})$ which is already
small and nonzero. For the last, let
\begin{align*}
a
&=
r_4 + r_5\\
&= -(r_1+r_2+r_3)\\
&= -z^{-4} - \frac13 z^{-2} + \frac1{81} z^2 - \frac{82}{243} z^4 + \frac43 z^6 - \frac{5585}{6561} z^8 + O(z^{10})
\end{align*}
[sorry, I switched from the earlier $a = -(r_4+r_5)$],
and proceed as we did before for $r_3$, finding
$$
3y^3 + (7a-5)y^2 + (3a^3-10a^2+2a+52)y + 10a^4-5a^3-22a^2-a-20
$$
with a somewhat more complicated power series
$$
82 z^6 - 34z^8 - \frac{58}{3} z^{10} - \frac{400}{3} z^{12} + O(z^{14}),
$$
but still enough for our purposes once $|y|$ exceeds $100$ or so.
The remaining candidates, $(X-r_1)(X-r_3)$ and $(X-r_2)(X-r_3)$, are tougher,
though related by $z \leftrightarrow -z$ so we can treat them together.
This time
$a = r_1+r_3$ is $z^{-4} + z^{-3} + \frac13 z^{-2} - \frac1{81} z^2 + \ldots$
involves powers of both $z^2$ and $z^3$, and it takes a while before
we can find a polynomial with no singular part. It's easier to use also
$b = r_1 r_3 = z^{-7} + \frac13 z^{-5} - \frac1{81} z^{-1} + \frac{82}{243} z
- z^2 + \ldots$,
which must also be an integer if $(X-r_1) (X-r_3)$ is to be an integral
factor of $q(X,y)$. Even then we need to use the monomials up to degree $22$,
$$
{}^{1,\; a,\; y,\; b,\; a^2,\; ay,\; ab,\; a^3,\; y^2,\; by,\; b^2,\;
a^2b,\; a^4,\; ay^2,\; aby,\; a^3y,\; y^3,\; ab^2,\; a^3b,\; by^2,\;
a^5,\; a^2y^2,\; b^2y,\; a^2by,\; b^3,\; a^4y,\; ay^3,\; a^2b^2}
$$
because the first few combinations we find, such as
$y^3 + by^2 + 2aby + 2ab^2 + (-a^3+a+2)b$, are identically zero.
But the combination with coefficients
$$
{}^{120,\; 44,\; -83,\; -159,\; 103,\; 82,\; 144,\; -42,\; -5,\; -36,\;
-9,\; -27,\; 9,\; 24,\; 86,\; -19,\; 24,\; 13,\; 42,\; 0,\; -14,\; 1,\;
-12,\; -10,\; -8,\; 9,\; 0,\; -9}
$$
succeeds, giving
$$
72z + 42z^2 - 214z^3 + 370z^4 - 134z^5 + \frac{142}{3}z^6 - \frac{1702}{9}z^7
+ \frac{2674}{27} z^8 + \frac{1064}{3} z^9 + O(z^{10}).
$$
The resulting upper bound on $y = z^{-6}$ still seems impractically large: certainly
nowhere near the astronomical $10^{20000}$ that the OP obtained
from the literature, but still a bit beyond the OP's computation to $10^9$
(it's about $72^6$, which exceeds $10^{11}$). However, we don't need
to try all $y$ up to that bound, only those which make our polynomial
nearly an integer (positive or negative, because replacing $r_1$ by $r_2$
yields the same power series with $-z$ in place of $z$); and there are
only a few hundred such to compute. For each of them, check whether
the nearest integer $y$ yields a reducible $p(X,y)$, and you're finally done.
Best Answer
The intuition for this method of passing from a rational solution to an integral solution seems pretty simple to me: passing from a rational solution to a nearby integral point (not necessarily a solution) is passing to a point whose denominators are 1, so you can anticipate that when you intersect the line through your rational solution and the nearby integral point with whatever curve or surface contains your solutions, the second intersection point on that line will have denominators that have moved closer to 1. That is, connecting a rational solution with some integral point will spit out a new solution whose denominators are somewhere between the denominators of your solution and the denominators of the integral point you used to produce the line.
Of course intuition is one thing and checking the details is another: you choose the integral point nearby and the math has to work out to show the denominators really get smaller in the second solution you produce. For instance, this method of proving the 3-square theorem goes through without a problem for a similar 2-square theorem (if an integer is a sum of two rational squares than it's a sum of two integral squares by the same method, replacing the sphere x^2 + y^2 + z^2 = a with the circle x^2 + y^2 = a). But this intuitive way of creating an integral solution from a rational solution breaks down if you apply it to the 4-square theorem: the inequalities in the proof just barely fail to work (sort of like doing division with remainder and finding the remainder is as big as the divisor instead of smaller).
The intuition also breaks down if you slightly change the expression x^2 + y^2 (sticking to two variables). Consider x^2 + 82y^2 = 2 and the rational solution (4/7,1/7). Its nearest integral point in the plane is (1,0), and the line through these intersects the ellipse in (16/13,-1/13), so the denominator has gone up. There actually are no integral solutions to x^2 + 82y^2 = 2. Or if we take x^3 + y^3 = 13 and the rational solution (2/3,7/3), its nearest integral point in the plane is (1,2), the line through these meets the curve again in (7/3,2/3), whose nearest integral point in the plane is (2,1), the line through them meets the curve in (2/3,7/3),...
A few years ago when I was giving some lectures on the method of descent, I worked out some examples of this geometric "three-square" theorem (start with an equation a = x^2 + y^2 + z^2 where a is an integer and x, y, and z are rational and produce in a few steps an equation where x, y, and z are integral) and I noticed in my initial examples that the denominators in each new step did not merely drop, but dropped as factors, e.g., if the common denominator at first was 15 then at the next step it was 5 and then 1. Maybe the denominators always decreas through factors like this? Nope, eventually I found a case where they don't: if you start with
13 = (18/11)^2 + (15/11)^2 + (32/11)^2
then the integral point nearest (18/11,15/11,32/11) is (2,1,3) and the line through these two points meets the sphere 13 = x^2 + y^2 + z^2 in the new point (2/3,7/3,8/3), so the denominator has fallen from 11 to 3, which is not a factor. (At the next step you will terminate in the integral solution (0,3,2).)