I would suggest you use the GWAS by GBLUP approach to make the problem a 500x500 size problem instead of 30K x 30K, where you can still recover the 30K effects with a simple back transformation. That model should take only few seconds to run instead of hours or days. See the last section of the following vignette:

https://cran.r-project.org/web/packages/lme4breeding/vignettes/lmebreed.qg.html