Having only a single missing value slows lmFit down by over an order of magnitude:

library("limma")
library("tictoc")

n_genes <- 10^6

sd <- 0.3*sqrt(4/rchisq(n_genes,df=4))
y <- matrix(rnorm(n_genes*6,sd=sd),n_genes,6)
y[1:2,4:6] <- y[1:2,4:6] + 2
design <- cbind(Grp1=1,Grp2vs1=c(0,0,0,1,1,1))

y_NA <- y
y_NA[1,1] <- NA

tic()
fit <- lmFit(y,design)
toc()

tic()
fit <- lmFit(y_NA,design)
toc()

While the first fit takes about 1.1sec, the second needs over a minute. Is this a bug?

The timings you give show that lmFit is actually very fast, especially so when there are no NAs. You are the first person ever to view that as a "bug".

lmFit does an intial scan for NAs or weights and, if they are absent, then it runs a special super-fast algorithm that only works when there are no NAs.