I want to analyse the relationship between a continuous variable "age" with the methylation of CpG sites.

To do this I created a design matrix to analyse in limma with design <- model.matrix(~ age, df)

But now I want to include a confounder such as sex (used as a factor). When I use design <- model.matrix(~ age + sex, df) I am unclear how to interpret the resulting topTable of significant sites as there is a columns of the table are "age, sexMale, AveExpr, F, P.Value, adj.P.Val".

What does the adj.P.Val indicate here? I have read the limma workflow and have not been able to determine the answer.

Greatly appreciate any help with this.

adj.P.Val is the p-value adjusted for multiple testing. If topTable() is run with default settings, then the adj.P.Val is the false discovery rate (FDR).

If you want to test for age adjusted for confounders,, then you need to run topTable with coef="age". From what you say, it sounds like you have left the coef argument blank. If you specify coef explicitly, then the resulting topTable will have exactly the same columns regardless of the number of confounders.