You have a highly oscillatory integrand — your integrand is oscillating around 30000 times in your integration domain, and you should expect that any generic quadrature scheme will need several quadrature points per oscillation, so you should figure on 100k+ quadrature points for just a 1d integral.
Indeed, calling QuadGK.jl explicitly to get an evaluation count, I find:
julia> using QuadGK
julia> quadgk_count(x -> integrand(x, 0), 0, 4500, rtol=1e-5, atol=1e-5)
(0.001980687766630197, 9.987036693025001e-6, 968415)
or almost 10^6 evaluation points for your 1d integral. There are ways to speed this up (e.g. increasing the quadrature order, pre-allocating buffers, etc.) but fundamentally this is a hard problem for a generic quadrature routine.
If your multidimensional integrand is similar along each dimension, you can see why — in 2d it would take \approx 10^{12} function evaluations.
You really need to re-think what you are doing — this is not something you are going to be able to brute force without exploiting more structure of your problem.
In certain cases, highly oscillatory integrals can be computed very efficiently — for example, if the oscillatory part is known and simple (e.g. a cosine), but it is multiplied by some complicated function (as in a Fourier transform), there are efficient methods. See e.g. [RFC/ANN] OscillatoryIntegralsODE.jl: Levin method + OrdinaryDiffEq and Integrate High Dimension highly oscillating integral