Hello,
As of v1.4.0, fmriprep computes the WM and CSF regressors from the BOLD pre-AROMA time series. In addition to the discussion in the issue linked in the OP, you can find simulation results here that evaluate potential implications of this implementation. There’s also a previous post/thread that briefly discusses this here.
I haven’t used denoiser, so I can’t make software-specific comments on usage, but removing the WM/CSF signal and filtering the data sounds reasonable. (We’ve also found that GSR can enhance ICA-AROMA performance.) In general, I would strongly encourage using a simultaneous bandpass+regression approach if at all possible (for instance, fitting a cosine basis filter to the data during the confound regression step) instead of regressing and then filtering.