Using ComBat-Seq to adjust for batch effect/ colony ID

[mastoldt]

I am currently working on an RNA-Seq dataset with an unbalanced design and two factors: treatment and colony ID (I am working with social insects, so colony ID is important to control for because of pseudoreplication).
Unfortunately, the samples had been dissected and extracted on three different days and those are in line with the colony ID.

When I plot a PCA, I clearly see an influence of colony ID/ extraction day, so a potential batch effect I would like to remove, although I cannot be certain whether this is due to extraction day or colony ID.

That's where I found out about ComBat-Seq and I thought this might solve the problem I have by removing the strong batch effect from my data. But I am unsure whether afterwards, I still need to include colony ID into my DESeq2 model, because I have already removed the batch effect which is the same as colony ID. But in my case individuals from the same colony represent pseudoreplicates, so not controlling for colony ID seems a bit odd to me.

So, should I run

a) ComBat-Seq + DESeq2 (gene expression ~ treatment) or
b) ComBat-Seq + DESeq2 (gene expression ~ colony ID + treatment)

Thank you a lot in advance, I am happy for any feedback on this.