Realize that the quadrature order N is not the total number of quadrature points — the number of quadrature points used for each subinterval is 2N+1, but the whole point of an adaptive algorithm like quadgk is that it subdivides the interval into more and more subintervals until convergence is achieved.
You want a relatively low-order rule for each subinterval so that it doesn’t waste integrand evaluations unnecessarily. In most cases you shouldn’t change the quadrature order from the default unless you are trying to use BigFloat precision with a smooth function and need to increase the convergence rate to get a huge number of digits.
It looks like the algorithm QuadGK.kronrod is using to get the Gauss–Kronrod points and weights is hitting the limits of Float64 precision if you ask it for the Gauss–Kronrod rule of order 1000, but normally one should never want this so I guess the algorithm was not designed with that kind of scaling in mind. (The same algorithm works fine for BigFloat precision, i.e. if you do QuadGK.kronrod(BigFloat, 1000) it succeeds.)
(If you need Gaussian-quadrature rules of extremely high order, see also the FastGaussQuadrature.jl package, though it’s not designed for adaptive integration, i.e. it doesn’t do Gauss–Kronrod or Gauss–Patterson weights.)