bf(A~ C + D+ scale(B) + (1|RE), family =negbinomial(link = "log")) +
bf(B ~ C+ D+ (1|RE), family =negbinomial(link = "log"))

The problem is that this model gives me an overall estimate of C and D on A and does not allow me to disentangle the direct and indirect effects. Is there a way to obtain estimates for this with brms or another package?

Have you tried to assess the same model using lavaan R package? And if you are looking for Bayesian version then you have blavaan Rpackage. The difference of course is that blavaan is based on Gibbs sampling and not Stan. The key advantage (if you are a R user) is that lavaan and blavaan give you the FIT indices, Modification Scores so you can assess your SEM systematically in a single go (which I would advice is the right way to do it).

If you want to use Stan/brms for SEM, I think you should approach in a Path Analysis approach where you build individual regressions along the paths for direct and indirect paths. However, note that you are not going to get any of the required SEM type output such as Fit Indices and Modification Scores.