Why is my quasibinomial GLM estimator biased - Monte Carlo simulation

Question

I'm playing with some Monte Carlo simulations to get an idea of the properties of some linear and non-linear models. The linear OLS model in my case is specified as:

$Y_t = \beta_0 + \beta_1x+ \epsilon,$

with $\beta_0=\beta_1=0.5$

My point of interest here is to figure out if the estimates $\hat{\beta_0}$ and $\hat{\beta_1}$ are unbiased and roughly equal to the true parameters $\beta_0$ and $\beta_1$.

Below, you can see my R code for doing a Monte Carlo simulation with both the linear OLS model and a non-linear quasibinomial probit GLM version of it.

#Monte Carlo simulation

N <- 100 #number of observations
mean <- 0 #mean
sd <- 1 #standard deviation
x <- rnorm(N,mean = mean,sd=sd) #simulate random exogeneous variables

n_sim <- 1000 #number of simulations

b0 <- 0.5 #true parameter beta 0
b1 <- 0.5 #true parameter beta 1

b0_lin <- numeric(n_sim) #stores the estimates of b0
b1_lin <- numeric(n_sim) #stores the estimates of b1
b0_probit <- numeric(n_sim) #stores the estimates of b0
b1_probit <- numeric(n_sim) #stores the estimates of b1

#apply the Monte Carlo simulation
for (i in 1:n_sim){
  set.seed(1+i)
  #simulate dependent variable
  y_lin <- b0+b1*x+rnorm(N, mean=mean, sd=sd) 
  #linear model
  mod <- lm(y_lin~x)
  #save vector of estimated parameters
  b0_lin[i] <- coef(mod)[[1]]
  b1_lin[i] <- coef(mod)[[2]]
  
  #transform dependend variable with inverse probit link funktion, thus it lies between zero and 1
  y_prob <- pnorm(y_lin)
  #quasibinomial probit GLM model
  mod <- glm(y_prob~x, family = quasibinomial("probit"))
  #save vector of estimated parameters
  b0_probit[i] <- coef(mod)[[1]]
  b1_probit[i] <- coef(mod)[[2]]
}

#Save means of the parameter vectors
b0_lin_mean <- mean(b0_lin)
b1_lin_mean <- mean(b1_lin)

b0_probit_mean <- mean(b0_probit)
b1_probit_mean <- mean(b1_probit)


Montecarlo_comparison <- data.frame("b0_mean"=c("True parameters"=b0,
                                               "Linear"=b0_lin_mean,
                                               "Probit"=b0_probit_mean),
                                   "b1_mean"=c(b1,
                                               b1_lin_mean,
                                               b1_probit_mean)
                                   )
print(Montecarlo_comparison)

As you can see below, the estimator from a linear OLS model is not biased. However, the quasibinomial probit GLM estimator appears to be severely biased.

My questions: Why is the quasibinomial probit GLM estimator biased in my Monte Carlo simulation? Is the process of my Monte Carlo simulation technically correct or did I maybe make a mistake in the model specifications? If I did it correct, what is the explanation for the fact that the quasibinomial probit GLM estimator is biased here?