Principal Value in NIntegrate. Integrand with multiple poles

I want to evaluate the numerical integral:

I1[y_] := NIntegrate[PolyGamma[(3*((0.61^2) - x (1 - x)*(y)^2))/(4*0.001) + 1], {x, 0, 1}]

The plot of the function defined by the integral above is:

The result is highly oscillatory for y>2*0.61. This happens because the PolyGamma function has poles if the argument is a negative integer or zero. My idea is to use Principal Value and exclude the points which the argument is negative and integer.

I reckon that the best approach is:

I2[y_] := NIntegrate[PolyGamma[(3*((0.61^2) - x (1 - x)*(y)^2))/(4*0.001) + 1], {x, 0, 1}, Method ->PrincipalValue, Exclusions -> {(3*((0.61^2) - x (1 - x)*(y)^2))/(4*0.001) + 1 == n}]

But I don’t know how to define "n" as a negative integer. I have two questions:

1. How can I define n as a negative integer and use it to exclude the problematic points?

2. There is another way to compute the principal value of this integral?

Integrate[(z/(k (k + z)) - 1/z - EulerGamma) /. 
 z -> Rationalize[(3*((0.61^2) - x (1 - x)*(y)^2))/(4*0.001) + 1, 
0], {x, 0, 1}, Assumptions -> y >= 0]

 (*-EulerGamma + 1/k - (
4 ArcTan[(50 Sqrt[3] y)/Sqrt[11203 + 40 k - 7500 y^2]])/(
5 Sqrt[3] y Sqrt[11203 + 40 k - 7500 y^2]) - (
4 ArcTan[(50 y)/Sqrt[11203/3 - 2500 y^2]])/(
5 y Sqrt[33609 - 22500 y^2])*)

F[y_, M_] := -EulerGamma + 
Sum[1/k - (4 ArcTan[(50 Sqrt[3] y)/Sqrt[11203 + 40 k - 7500 y^2]])/(5 Sqrt[3] y 
Sqrt[11203 + 40 k - 7500 y^2]), {k, 1, M}] 
- (4 ArcTan[(50 y)/Sqrt[11203/3 - 2500 y^2]])/(5 y Sqrt[33609 - 22500 y^2])

Plot[F[y, 2000] // Re, {y, 0, 2}, PlotRange -> {Automatic, {2, 6}}]

Plot[F[y, 5000] // Re, {y, 0, 2}, 
PlotRange -> {{1.225, 1.5}, {2, 6}}](*In range 1.225 < x < 1.5*)