> my_model <- gls(variable~Soil_type+Genotype+Soil_type:Genotype, data=my_data, na.action = na.omit)

> shapiro.test(resid(my_model, type = "normalized"))

    Shapiro-Wilk normality test

data:  resid(my_model, type = "normalized")
W = 0.99199, p-value = 0.8568

> bartlett.test(resid(my_model, type = "normalized") ~ fitted(my_model, type = "normalized"))

    Bartlett test of homogeneity of variances

data:  resid(my_model, type = "normalized") by fitted(my_model, type = "normalized")
Bartlett's K-squared = 29.118, df = 29, p-value = 0.4589

> anova(my_model)
Denom. DF: 62
                   numDF   F-value p-value
(Intercept)            1 1664.3700  <.0001
Soil_type              2  121.1435  <.0001
Genotype               9    3.9401  0.0005
Soil_type:Genotype    18    1.3449  0.1930