As in @MichaelLew's suggestion, the two samples
of differences (After - Before) for each group might be similar
to the fictitious data sampled in R below. (In the plots of your question, it seems differences might be roughly
normal, but with a larger sample SD for the second
sample)
set.seed(2022)
d.1 = rnorm(500, 0, 1.5)
summary(d.1); sd(d.1)
summary(d.1); sd(d.1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-5.21285 -1.10637 -0.05737 -0.04802 0.93023 5.64779
[1] 1.527096 # sample SD
d.2 = rnorm(550, 3, 2.5)
summary(d.2); sd(d.2)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.966 1.394 2.906 3.119 4.953 10.682
[1] 2.488868
boxplot(list(d.1,d.2), horizontal=T,
col=c("skyblue2","wheat"))

Because the two samples of differences may not
have the same variances, a Welch 2-sample t test
seems appropriate. The P-value for a 2-sided test
is nearly $0,$ so the difference between the two
samples is significant at any reasonable level.
t.test(d.1, d.2)
Welch Two Sample t-test
data: d.1 and d.2
t = -25.096, df = 923.59, p-value < 2.2e-16
alternative hypothesis:
true difference in means is not equal to 0
95 percent confidence interval:
-3.414837 -2.919489
sample estimates:
mean of x mean of y
-0.04802171 3.11914101
If you doubt that your real data are nearly normal,
then you could use a nonparametric Wilcoxon rank
sum test---instead of a Welch t test.
However, the if the two samples are of different shapes (including different variations, as for my fictitious data),
the Wilcoxon test should not be regarded as a test
of different medians, but as a test of stochastic
domination. Again the difference between the locations of the two samples is highly significant.
wilcox.test(d.1, d.2)
Wilcoxon rank sum test
with continuity correction
data: d.1 and d.2
W = 37731, p-value < 2.2e-16
alternative hypothesis:
true location shift is not equal to 0
Roughly speaking, d.2
tends to have larger
values than d.1
. This can be shown by comparing
the empirical CDFs (ECDFs) of the two samples: points for d.2
(brown) tend to be to the right (and thus
below) points for d.1
.
plot(ecdf(d.2), main="ECDF Plots", col="brown")
lines(ecdf(d.1), col="blue")
