Getting different p values for main effects alone than when in a complex model

[mwhite942]

I am trying to use adonis2 in RStudio to run a PERMANOVA on a Bray-Curtis matrix for microbiome sequencing data (OTUs). I have multiple independent variables of interest, and so I would like to run a complex model: model1 <- adonis2(matrix ~ rep + line*loc*compartment, data = metadata1, permutations = blocking) where "blocking" is defined as blocking <- how(within = Within(type = 'free'), plots = Plots(type = 'none'), blocks = metadata1$rep, nperm = 9999, observed = T) in order to account for my randomized complete block design. I get the following results:
`Terms added sequentially (first to last)
Blocks: metadata1$rep
Permutation: free
Number of permutations: 9999

adonis2(formula = matrix ~ rep + loc * compartment * line, data = metadata1, permutations = blocking
Df SumOfSqs R2 F Pr(>F)
rep 7 5.47 0.01482 51.5374 0.0001
loc 3 197.03 0.53362 618.7016 0.0001
compartment 2 18.81 0.05095 88.6026 0.0001
line 11 2.76 0.0748 2.3646 0.0001
loc:compartment 6 30.04 0.08135 47.1597 0.0001
loc:line 33 5.06 0.01369 1.4434 0.0001
compartment:line 22 1.95 0.00527 0.8339 0.9382
loc:compartment:line 66 4.83 0.01308 0.6896 1.0000
Residual 973 103.29 0.27973
Total 1117 369.24 1.00000`

When I make NMDS plots, I can clearly see that the points separate by loc (location), and to a lesser extent compartment and rep, so significant results for those make sense. I am a bit surprised that compartment and rep have p values that are just as significant as the loc p value, but I understand that this may be a quirk of the PERMANOVA (ie, maybe I would need to increase the permutations a lot to get different p values). What I am more concerned about is the result of the "line" variable. In the NMDS plots, there is no visible clustering by line, and when I add 95% CI ellipses, the ellipses, while not being exactly on top of each other, a ll encircle nearly all of the points and very very nearly completely overlap. When I run an adonis2 with just the line variable (matrix ~ line), I get a p value of 0.94 with 99 reps (I only did a low number of reps because I wanted to shorten the run time; similar tests with only the other variables did give significant p values). So, I am wondering why I am be getting such vastly different p values when I include the "line" variable in the complex model than when I run it on it's own. Is this a coding problem, or some kind of bug? Am I not understanding how the test works correctly? I am also (less) concerned that, when I run the full model, all of my variables except for the last 2 interactions always give identical p values to each other, despite having wildly different F values. (I have run the model more than once, with the same result.) In case it's helpful/relevant, I also ran a betadisper test on each of the variables (I couldn't figure out how to run one on the entire model), and found that loc and compartment have significantly different variations within their groups, but line and rep both returned a p value > 0.9.