  In this work, we investigate a quasilinear subdiffusion model which involves
a fractional derivative of order $\alpha \in (0,1)$ in time and a nonlinear
diffusion coefficient. First, using smoothing properties of solution operators
for linear subdiffusion and a perturbation argument, we prove several
pointwise-in-time regularity estimates that are useful for numerical analysis.
Then we develop a high-order time stepping scheme for solving quasilinear
subdiffusion, based on convolution quadrature generated by second-order
backward differentiation formula with correction at the first step. Further, we
establish that the convergence order of the scheme is
$O(\tau^{1+\alpha-\epsilon})$ without imposing any additional assumption on the
regularity of the solution. The analysis relies on refined Sobolev regularity
of the nonlinear perturbation remainder and smoothing properties of discrete
solution operators. Several numerical experiments in two space dimensions show
the sharpness of the error estimate.
