Can metagenomSeq’s fitZig calculate coefficient for interaction terms (eg; gender*birthweight)?

Yes, just when you are defining your model.matrix, explicitly include the terms.

i.e.

mod = model.matrix(~gender*birthweight,data=phenodata)