How to describe the relationship between these two variables? I would like to quantitatively describe the relationship between the two variables shown in this plot, but I am not sure what would be the correct way to do so. More specifically, my aim would be to assess whether there is a relationship between both variables and to quantify the strength of such a relationship (e.g. quantifying the amount of V2 variance explained by V1 with the coefficient of determination R2 or any other suitable statistic)

I have tentatively used the GAM model available in the geom_smooth function of ggplot  (method="gam",formula = y ~ s(x, bs = "cs")
and I have reproduced the same model with the gam package (see code and output below). However, I am not sure if this would be enough to claim that there is a relevant relationship between V1 and V2 and, if so, what would be the crucial statistics that I should report to substantiate such a possible conclusion. On the other hand,  I have calculated the Spearman correlation trying to assess the relationship at the ordinal level (rho=0.480, p= 7.221e-08).
Thanks in advance for any help/ advice

library(gam)
v=gam(V1~s(V2),bs="cs", data=ff)

Call: gam(formula = V1 ~ s(V2), data = ff, bs = "cs")
Deviance Residuals:
     Min       1Q   Median       3Q      Max 
-3.36459 -0.31432 -0.14645  0.09815  7.18685 

(Dispersion Parameter for gaussian family taken to be 1.7516)

    Null Deviance: 377.6946 on 115 degrees of freedom
Residual Deviance: 194.4251 on 111.0001 degrees of freedom
AIC: 401.1024 

Number of Local Scoring Iterations: 2 

Anova for Parametric Effects
           Df  Sum Sq Mean Sq F value    Pr(>F)    
s(V2)       1  71.352  71.352  40.736 4.159e-09 ***
Residuals 111 194.425   1.752                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova for Nonparametric Effects
            Npar Df Npar F     Pr(F)    
(Intercept)                             
s(V2)             3 21.299 5.709e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#############
############
EDIT:
Following the advices received  I have transformed V1 using log (both proposed transformations produced virtually identical results). Then, I have estimated the relationship through linear regression. However, the obtained fit does not seem to be much better...

 A: When bivariate data present this shape — here, with the points falling nearly vertically around x = 0 and x = 1, and then nearly horizontally when x > 2 — a continuous model is unlikely to produce satisfactory results.
Sometimes a segmented model, such as a linear-plateau or quadratic plateau, will work well enough.
There is, however, another approach called a Cate–Nelson model which may be helpful in these cases.  A Cate–Nelson approach simply tries to divide the data into a group with low x values and low y values, and into a group with high x values and high y values.†  Data that fall into these groups are considered to conform to the model.
You can find this approach used in agronomic studies.  There are different ways to determine the x value and y value used to divide these groups. One could simply use an iterative approach to find the values that best divide the data, or one could just estimate the values visually. Cate and Nelson (1971)†† propose a method that divides the data to maximize the sum of squares.
I included a plot that roughly mimics the data in the post.  Here, the data are divided at c. x = 1.15 and y = 0.62, and 88% of points fall into quadrants II and IV.

† Of course, the model could just as well group data into low x values with high y values, and so on.
†† Cate, R. B., & Nelson, L.A. (1971). A simple statistical procedure for partitioning soil test correlation data into two classes. Soil Science Society of America Proceedings 35, 658–660.

A: The first thing that stands out is that V1 looks highly right skewed. A first pass would be to log-transform V1 and see if a linear regression then provides a reasonable fit.
Given that the data contain zeros, a common alternative to log transformation that preserves zeros is the inverse hyperbolic sine:
ihs = log(V1 + sqrt(V1^2 + 1))

See, e.g.:
Burbidge, John B., Lonnie Magee and A. Leslie Robb. 1988 "Alternative
Transformations to Handle Extreme Values of the Dependent Variable."
Journal of the American Statistical Association, vol. 83, 123-127.
A: If you want "to assess whether there is a relationship between both variables and to quantify the strength of such a relationship," then any general rank-correlation coefficient, for example Spearman's or Kendall's, would be a reasonable choice. That makes no assumptions about the underlying distributions, the functional form of the relationship, or which of your 2 variables is a predictor versus an outcome. You have a highly significant Spearman coefficient, consistent with a strong relationship between them.
The Spearman coefficient helps with "quantifying the amount of V2 variance explained by V1," if you are willing to consider the variance in ranks rather than variances in the actual values: the Spearman coefficient is the correlation coefficient between the 2 sets of ranks.
For "quantifying the amount of V2 variance explained by V1" in their original scales, then regression is an approach. Given the non-linear relationship between them, however, you presumably want to choose some transformation(s).
The best transformation(s) depends on how you intend to use the result. If you don't have a specific functional form in mind based on your understanding of the subject matter, generalized additive models (GAM) as you tried (although you modeled V1 as a function of V2, not the direction you seem to be interested in) and spline-based approaches provide general ways to proceed. Note, however, that you might tend to overestimate "the amount of V2 variance explained by V1" with those approaches: the results might not generalize well to new data samples unless you took precautions like penalization to minimize the effects of vagaries in fitting the particular data sample at hand.
If you don't care about the functional form of the relationship between the 2 variables, then the Spearman correlation is a very good choice. You can think of it as their correlation following one of the simplest of all transformations: each is simply transformed into its empirical cumulative distribution function.
